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Abstract 



In the study of condensed matter systems, a paradigm which is very generally invoked is that of 
a ground state and its low energy excitations. A classical solid, for example, is a periodic crystal 
with low energy, long wavelength phonon fluctuations. Often these low energy excitations appear 
due to the complete or partial breaking of a spatial symmetry e.g. rotational or translational 
symmetry. What happens if explicit constraints are introduced such that these low energy modes 
are unavailable? 

This question has assumed some importance in recent years due to the advent of nano technology 
and the growing use of nanometer scale devices and structures. In a small system, the size limits 
the scale of the fluctuations and makes it imperative for us to understand how the response of 
the system is altered in such a situation. In this thesis, this question is answered for the special 
case of interfacial fluctuations in two dimensions (2d). The energy of an interface between two 
phases in equilibrium is invariant with respect to translations perpendicular to the plane (or line 
in 2d) of the interface. We study the consequence of breaking this symmetry explicity using an 
external field gradient. One expects that since low energy excitations are suppressed, the interface 
would be flat and inert at all times. We show that surprisingly there are novel fluctuations and 
phenomena associated with such constrained interfaces which have static as well as dynamic 
consequences. 

After a couple of chapters containing introductory material on interfaces, in Chapter 3 we inves- 
tigate static and dynamic properties of an Ising interface defined on a 2d oriented square lattice. 
The interface is constrained by a non-uniform external field with a profile which is such that 
the field changes sign across a (straight) line which we call the "edge" of the profile. The Ising 
interface lies as close to the edge as possible and the presence of the symmetry breaking field 
suppresses long wavelength interfacial fluctuations as expected. The interface therefore remains 
essentially flat over long length scales. At short scales of the order of a few lattice parameters, 
however, the interface begins to show "self-assembled" patterns depending on the orientation 
of the edge. These patterns may be indexed by all possible rational fractions. The energy of 
the interface depends on the index of the rational fraction and the interface deforms locally to 
conform to low energy states. The local orientation of the interface plotted as a function of 
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Abstract 



the average orientation shows a complete Devil's staircase structure. If the external field profile 
is now translated perpendicular to the interface with a fixed velocity, the interface follows the 
profile, but we now have transitions among the patterns which occur in a way so as to render 
the Devil's staircase incomplete. Low index patterns are stable while patterns indexed by higher 
order rationals either transform to low index ones or remain "amorphous" i.e. fluctuate randomly 
between nearby stable states. At very high velocities the interface detaches from the field edge 
and coarsens dynamically. We believe that our results have some relevance on real crystal growing 
techniques apart from being of theoretical interest due to the presence of an infinite hierarchy of 
driven dynamical transitions. 

The later part of the thesis deals with a more realistic interface (also in 2d) namely that between a 
classical liquid and a solid. Similar to the Ising interface, this is created using a chemical potential 
field in the form of an atomic trap. The presence of the trap allows the atoms to condense into 
a solid which is separated from the surrounding liquid by a liquid solid interface. For reasons 
similar to the Ising case, the interface is flat and "crystallization waves" which normally tend to 
roughen a solid liquid interface are suppressed. The small size solid thus formed can exchange 
particles with the liquid through the interface. We show that this transfer of particles happens 
predominantly by a coherent addition and removal of entire atomic layers as the depth of the 
potential well is changed. The elastic energy cost of flaking off incomplete layers is prohibitively 
large if the thickness of the crystal is below a certain critical amount. This ensures that only 
whole, single, layers are exchanged. This provides a novel mechanism of stress relaxation in a 
nano-sized single crystal. Topological defects like dislocations tend to destroy this coherence but 
are rapidly eliminated once they appear on the interface as steps. 

Finally, we study momentum and energy transfer across the liquid solid interface in the presence of 
this "layering" transition using molecular dynamics simulations. With the depth of the potential 
held close to the value required for a layering transition, a tiny momentum pulse is generated, 
near a free surface of the crystal. The pulse travels as a shock wave through the solid and 
emerges on the other side taking with it an entire crystal layer; the ejected layer carrying away a 
definite amount of momentum and energy behaving, as it were, like a coherent "quasi particle" 
in a completely classical context. This results in a prominent peak in the absorption spectrum. 
We believe that this driven layer ejection may have practical applications apart from its interest 
from a fundamental viewpoint. The transfer of energy through the interface is characterized by 
the Kapitza or contact resistance which we measure using non-equilibrium molecular dynamics 
where a temperature gradient is setup using proper boundary conditions. We observe that the 
Kapitza resistance shows a discontinuity at the layering transition. This may enable one in future 
to design interfaces with controllable thermal properties 
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1 Introduction to interfaces and 
summary of the thesis 



To a large extent, the properties of most everyday materials are determined by the properties of 
their internal interfaces |2J- This is equally true, for example for a block of steel, where grain 
boundaries interact with dislocations to produce hardness and resistance to plastic deformation, or 
a catalytic convertor in an automobile exhaust where interface properties control crucial chemical 
reactions, or in a mammalian lung tissue where lung surfactants residing in the inner lining of 
alveoli reduce the amount of work required in breathing and thereby contribute directly to the 
survival of the organism. 



It is therefore obvious that any means that we discover, which alters the properties of interfaces 
will have a profound effect on the overall structure and properties jof many materials of tech- 

which attempts precisely 



nological importance. Indeed, a large amount of work exists j^, Q, 
that. Interfacial properties may be altered using alloying elements |6 
prevents segregation of specific elements at interfaces, surfactants 
tension, stress etc. 



which either segregate or 
which reduce the surface 



In this thesis, we explore the possibility of tuning interfacial structure and properties using 
non-uniform external fields. We show that large potential gradients at interfaces reduces the 
dominant low-energy, long wavelength fluctuations. Instead of such interfaces becoming inert, 
however we show that novel fluctuations are generated which lead to interesting behaviour which 
may have some technological applications. Such fine tuning of interfacial fluctuations may be of 
use in the emerging area of nanotechnology j3| where the properties of interfaces at the length 
and time scales of interest in this thesis, is crucial. 

We begin this chapter by providing a brief introduction to the thermodynamics of interfaces 
followed by a discussion of the dominant low energy equilibrium fluctuations. We then introduce 
the basic paradigms for understanding the dynamics of driven interfaces. We end this chapter 
with a systematic overview and summary of the rest of the chapters of this thesis. 
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Introduction to interfaces and summary of the thesis 



1.1 Definition of interfaces and interfacial region 

By definition, bulk phases are regions of configuration space associated with a global minimum 
of an appropriate free energy. Often, instead of a single phase, one obtains a mixture of two or 
more phases which coexist. The region in-between two coexisting bulk phases is an interface, a 
surface being simply an interface where one of the bulk phases happens to be a gas or vacuum. 
An interface can also be composed of a material that is different from the two bulk phases such 
as surface coatings, epitaxial layers and membranes composed of amphiphilic molecules that may 
lie between bulk domains of water and oil. Although one may think of an interface as being 
of negligible thickness, in fact the thickness of the interfacial region is significant and definitely 
non-zero, when we consider phenomena occurring at nano-meter scales. 

Any intensive property, such as density, varying from, say, its value in a liquid phase through 
the interface to a gas phase gives a plot such as that in Fig. 1.1. In this case the density shows 




height 



Figure 1.1: Variation of density from the high density liquid to the lower density gas phase. Dashed 
lines indicate the interfacial region 

a smooth transition from the high density of the liquid to the much lower density of the gas. 
The bulk phases can be separated from the interface by two surfaces parallel to one another and 
positioned so that the bulk phases are homogeneous and uniform (uniform density in this case) 
while the inhomogeneity and non-uniformity are connected entirely within the interfacial region 
lying between the two surfaces. The dashed lines in Fig. 1.1 illustrate this point. However we 
should note that for some properties the transition from one bulk phase to another does not 
follow a smooth monotonic form such as that in Fig. 1.1. Instead we may have a behavior as 
shown in Fig. 1.2. As an example, the concentration of some solutes at the gas-solution interface 
may reach values much higher than those in either bulk phase and exhibit a profile such as that 
in Fig. 1.2(a). Such quantities are said to show an interfacial excess. Conversely, there may be 
relative depletion of the concentration of a solute at an interface (Fig. 1.2(b)) 



1.2 The thermodynamics of interfaces 
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Figure 1.2: (a) Variation of an extensive thermodynamic quantity with distance across an interface 
between X and y phases. The layer between the two dashed lines includes the interfacial region. Also 
shown is the "dividing surface" of Gibbs (b) Variation of the same except that here there is a relative 
depletion as opposed to interface excess in (a). 

1.2 The thermodynamics of interfaces 

Given the above definition of an interface it is possible to define interfacial analogues of bulk 
thermodynamic quantities. Formally, this may be achieved in several equivalent ways. Below, we 
follow the method of Cahn j^. An alternative derivation was formulated by Gibbs more than a 
century ago jlf. 



1.2.1 The interface free energy : Cahn's method 

Consider a system consisting of two bulk phases, X and 3^, each containing C components, which 
are in contact along a flat interface. The interface is growing in equilibrium within a container 
in contact with suitable reservoirs. The entire system, including the reservoirs, is maintained 
at constant temperature, T, hydrostatic pressure, V, and chemical potential, yUj, for each of 
the components. The increase in internal energy, d£, due to the accretion of dj\fi particles of 
component leading to corresponding changes in the entropy S, volume V and area A is given by, 

c 

dS = TdS -VdV + ^i2idAfi + afdA (1.1) 

1=1 

Eqn. 1.1 contains the usual bulk term plus the added term ajdA which is required to account 
for the increase in internal energy of the system which is associated with the increase in the area 
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of the interface. From Eqn. 1.1, we have 



d£ 



;i-2) 



and af, the interfacial free energy per unit area, is therefore defined as the increase in internal 
energy of the entire system per unit increase in interface area at constant S and V and Mi are 
held constant. We may as well define ct/ in terms of the Gibbs free energy, Q, the Helmholtz free 
energy, JF, and the grand potential, VL, as 







'dT' 




'dVL' 




M. 













;i.3) 



Further, on integrating Eqn 1.1, we obtain finally 



;i-4) 



where Q = £ + W — TS Since the quantity Yl^=i f^iM'i is the total Gibbs free energy which the 
homogeneous X and y phases would possess together, we may identify aj simply as the excess 
Gibbs free energy of the entire system per unit interface area due to the presence of the interface. 
Note that the entire contribution to af arises from the excess Gibbs free energy of the interfacial 
region. There is no contribution from the bulk. This ensures that aj is independent of our choice 
for the width of the interface as long as it is large enough to contain the full variation of all 
intensive interfacial quantities e.g., the region between the dashed lines in Fig. 1.1. 



1.3 Equilibrium fluctuations of interfaces 

Interfaces in equilibrium rarely remain flat. As can be easily seen, translating a planar interface 
between two coexisting phases in equilibrium in a direction perpendicular to the interface costs no 
energy. It therefore follows that long wavelength fluctuations of the interface would be energeti- 
cally cheap and would therefore dominate and roughen the regions between two materials. How 
strong these effects are depends on dimensionality (i.e., if the interface is a line or a plane) and 
on temperature. For instance, for fluid interfaces, two-dimensional systems with one-dimensional 
interfaces have a mean-square roughness of the interface that varies linearly with system size; 
the corresponding quantity in three-dimensional systems is logarithmic with the system size. To 
understand the role of fluctuations in interfaces, we consider first the simplest case where the fluc- 
tuations are thermal in origin. For fluid interfaces, in the limit of zero viscosity, these fluctuations 
are known as capillary waves. 



1.3 Equilibrium fluctuations of interfaces 
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1.3.1 Thermal fluctuations 

Consider the fluctuation 0, of a two-dimensional surface defined in the Monge representation 
as z = h{x,y). The area of the flat surface is denoted by A. For slowly varying fluctuations 
of this surface about a flat shape {h = h^, where is a constant) the additional surface free 
energy of the undulated interface over that of the flat one (AjF, = JF, — 7^) is approximately 

A^^ = ^7/^xrf,((|)^ + (|)^) (1.5) 

Note that the free energy AJF, is independent of h itself and depends only on the derivatives, 
due to translational invariance. In Fourier space, 

^^s = \lY.^'\h{c^)? (1.6) 
q 

The mean-square value of the fluctuating variable /i(q), using equipartition, can be written as 

<IMq)r>=^ (1.7) 

This gives the mean-square value of each Fourier mode in thermal equilibrium. 

The mean-squared, real-space fluctuations of the surface about a flat profile, are given by 

< h\r) >=\y< IMq)r >= tA^ /rfq < IMq)r > (1-8) 



A^ i-v^.i (2 



q 



To avoid logarithmic divergences (when q is two-dimensional) while performing the integral, we 
introduce a lower limit to the integral, vr/L, related to the finite size of the system (L oc \/~A) 
and an upper limit, tx / a due to the finite molecular size (proportional to a, the intermolecular 
separation). Then we have 

< h\v) >= -^log- (1.9) 
2777 a 

By symmetry, < /i^(r) > is independent of the point r in the x — y plane at which the fluctuation 
is calculated. Note that the mean-square fluctuation diverges as the logarithm of the system size. 
For a one-dimensional interface (i.e., a line instead of a surface) this divergence is even more 
severe and increases linearly with the size of the system. Thus, due to reduced dimensionality 
of these interfaces, the thermal fluctuations can have drastic effects on the "flatness" of the 
interfaces. 

1.3.2 Capillary waves 

For fluid interfaces, in the limit of negligible viscosity, the dominant fluctuations are capillary 
waves. We first write down the potential energy of an incompressible fluid surface acting under 
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gravity and surface tension. The free energy, JF,, is given by a constant term that comes from 
the bulk (where density variations are not allowed) free energy, Tq, and a contribution from the 
surface tension and gravitational terms : 

= + ^ y dA\p^gh^ + 7|V/ip] (1.10) 

where A is the area, po is the density, and 7 is the surface tension. Solving the equation of 
motion for the interfacial position h(x,y) resulting from the above free energy in the undamped, 
zero viscosity, inertial limit, we can obtain the dispersion relation Q 




u^=\q[9 + — (1.11) 



which shows that the frequency of such capillary fluctuations, ^ 0, as g ^ 0. For a liquid 
solid interface, these capillary fluctuations are known as "crystallization waves" jl^]. For most 
ordinary solids, such waves are highly damped, however, though they have been experimentally 
studied in great detail for liq helium at temperatures 0.03 to 0.5 K [12. Il3l| 



1.4 Morphology and dynamics of growth of interfaces 

in most situations, interfaces are actually not in equilibrium and it is difficult to define equilibrium 
structural and thermodynamic quantities. In such cases however, certain scaling laws, which 
determine how non-equilibrium interfacial fluctuations behave with time and system size come 
in useful in characterizing the structure of an interface. The scaling exponents behave similar to 
the exponents well known from studies of equilibrium critical point phenomena such that they 
can be used to define distinct "universality classes" of non-equilibrium, driven interfaces. In 
this section, we discuss briefly two such universality classes, viz. the Edwards-Wilkinson and the 
Kardar-Parisi-Zhang. 



1.4.1 Scaling concepts 

Consider a one-dimensional interface between two phases in two dimensions. Let us consider a 
situation where one of the phases is growing at the expense of the other. To describe the growth 
of this one-dimensional interface quantitatively, we introduce two functions (lO| |. 

• The mean height the surface, h, is defined by 

h{t) = y[ dxh{x,t) (1.12) 
L Jo 



1.4 Morphology and dynamics of growth of interfaces 
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where h{x,t) is the position of the interface at x at time t. The mean height increases 
linearly with time, 

h{t)r^t (1.13) 

if the interface moves at a constant rate. 

The interface width, which characterizes the rougliness of the interface, is defined by the 
rms fluctuation in the height. 



1 ''^ 



w{m)^\j- [hix,t)~hit)]^ (1.14) 

The time evolution of the surface width has two regions separated by a "crossover" time 

(i) Initially, the width increases as a power of time, 

w{L,t)r^tf^ [t«t^] (1.15) 

The exponent /3, which we call the growth exponent, characterizes the time-dependent dynamics 
of the roughening process. 

(ii) This is followed by a saturation regime during which the width reaches a saturation value, 
Wsat- As L increases, the saturation width, Wsat, increases as a power law, 

w,,i(L)~L- [t»Q (1.16) 

The exponent a, called the roughness exponent, characterizes the roughness of the saturated 
interface. 

(iii) The crossover time (or saturation time) at which the interface crosses over from the 
behavior of (i) to that of (ii) depends on the system size, 

tx ~ (1.17) 

where z is called the dynamic exponent. 

The exponents a, (3 and z are not independent. If we approach the crossover point (t^,w(t^)) 
from the left, we find, according to Eqn. 1.15, that w{tx) ~ tf . However, approaching the same 
point from the right, we have from Eqn. 1.16, w{tx) ~ L". From these two relations follows 

~ which according to Eqn. 1.17, implies that 

a 

. = - (1.18) 

Eqn 1.18, a scaling law linking the three exponents, is valid for any growth process that obeys 
the Family-Vicsek scaling relation |l^ 

w{L,t) Ly (1^ (1.19) 
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1.4.2 The Edwards-Wilkinson equation 

Consider an interface, defined by a height function h(x,t) which obeys the following stochastic 
growth equation, 

^^^ = G(/.,x,t) + r/(x,t) (1.20) 

where ?7(x, t) is a random noise and G{h, x, t) contains the deterministic part. Consider also that 
the equation of motion Eqn. 11.201 obeys the following symmetries |lo| : 

(i) Invariance under translation in time. The growth equation is invariant under the transfor- 
mation t ^ t + St- This symmetry rules out an explicit time dependence of G. 

(ii) Translation invariance along the growth direction. The growth rule is independent of where 
we define h = 0, so the growth equation should be invariant under the translation h ^ h + 5h- 
This symmetry rules out the explicit h dependence of G, so that the equation must be constructed 
from combinations of V/i, V^/i, V"/i, etc. 

(iii) Translation invariance in the direction perpendicular to the growth direction. The equation 
is independent of the actual value of x, having the symmetry x — > x + (5a;. This exclude explicit 
X dependence of G. 

(iv) Rotation and inversion symmetry about the growth direction n. This rules out odd order 
derivatives in the coordinates, excluding vectors such as V/i, V(V^/i), etc. 

(v) Up/down symmetry for h. The interface fluctuations are similar with respect to the mean 
interface height. This rules out even powers of h, terms such as (V/i)^, (V/i)"^, etc. A moments 
reflection tells us that this symmetry implies that h{x,t) represents in this case an interface 
between two equilibrium coexisting phases. 

The final form of the growth equation, consistent with all the symmetries listed above is, 

= iy^h) + iy% + ... + iy^^h) + iy^h)iyhf + ... 

at 

+{V^^h){Vhf^ + 7]{^,t), (1.21) 

where n,k,j can take any positive integer value. For simplicity of notation we do not explicitly 
indicate the coefficients in front of the terms. 

Since we are interested in the scaling properties, we focus on the long-time {t oo), long- 
distance (x — i> oo) behavior of functions that characterize the surface. In this hydrodynamic 
limit, higher order derivatives should be less important compared to the lowest order derivatives, 
as can be confirmed using scaling arguments. 

The noise term, 77(x, t) in Eqn. 1.21 incorporates the stochastic character of the fluctuation 
process and is assumed to be uncorrelated, with the properties 

(r?(x,t)) = (1.22) 
< r/(x,t)r7(x',t') > = 2D5'^(x - x')(5(t - t') (1.23) 



1.4 Morphology and dynamics of growth of interfaces 
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Thus, the simplest equation descri bing the fluctuations of an equilibrium interface, called the 
Edwards-Wilkinson (EW) equation [14] has the form 

^^^ = uV'h + v{^,t) (1.24) 

Here u is sometimes called a "surface tension", for the uV'^h term tends to smoothen the 
interface. Eqn 1.24 is valid in the small gradient approximation, i.e., in the limit (V/i) << 1. 
The average velocity of the interface is zero, which can be seen by using Eqn 1.24, 

The contribution from the Laplacian term is zero, due to the periodic boundary conditions, and 
since the noise has zero average according to Eqn 1.22. 

The scaling exponents of the EW equation can be obtained either by using scaling arguments 
or by solving the equation. The exponents in d-dimensions are given as follows : 

2-d ^ 2-d 
« = ^, /^ = ^' ^ = 2 (1-26) 

Thus for (i = 2, we find that a = P = 0, i.e., the width scales logarithmically with time at 
early times, and the saturation width depends on the logarithm of the system size. For d > 2, 
the roughness exponent a becomes negative, which means that the interface is flat, because the 
surface tension suppresses any noise-induced irregularity. 

1.4.3 Non-equilibrium interface : Kardar-Parisi-Zhang equation 

The first extension of the EW equation to include non-linear terms was proposed by Kardar, Parisi 
and Zhang (KPZ) [l^, The KPZ equation is the simplest growth equation which has the 
symmetries (i)-(iv) of the linear theory discussed in Eqn. 1.21, but the up-down symmetry of 
the interface height h(x,t) is broken. The source of the symmetry breaking is, of course, the 
existence of a driving force, F, perpendicular to the interface, and uniform in space and time 
which selects a particular growth direction for the interface. In the linear theory this up-down 
symmetry excludes terms such as (V/i)^" from the growth equation. The lowest order term of 
this sort is (V/i)^ which, if added to the EW Eqn. 1.26, results in the KPZ equation given as 

^^^ = uV'h+^{Vhr + r^{x,t) (1.27) 

In general, the existence of a driving force is a necessary, but not sufficient condition for the 
broken up-down symmetry in h, and hence for the appearance of the nonlinear term. As an 
example, consider the random deposition (RD) model with surface relaxation Jl^. In the RD 
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model, each particle falls along a single column toward the surface until it reaches the top of the 
interface, whereupon it sticks irreversibly. To include surface relaxation, the deposited particle is 
allowed to diffuse along the surface up to a finite distance, stopping when it finds the position 
with the lowest height. In this model, particle deposition generates a driving force that makes 
the interface grow. Deposition apparently breaks the up-down symmetry of the growth, but the 
model is nevertheless described by the EW equation, so in fact no symmetry breaking occurs. If 
we transform to a system of coordinates moving together with the interface, the model generates 
an interface which is invariant under the transformation h —h. 

Scaling exponents for the KPZ equation can be obtained exactly for d = 1 via a dynamical 
renormalization group calculation, to be a = 1/2, P = 1/3 and z = 3/2. However for d > 1, 
dynamic RG analysis does not give the scaling exponents. There are several models that predict 
the scaling exponents in higher dimensions but none of them are exact. For d = 2, large scale 
simulations predict that (3 = a/z = 0.240 ± 0.0001 (see jl^ and references therein). There is 
however a scaling relation that is true in any dimensions : 

a + z = 2 (1.28) 



1.5 Plan of the thesis 



aces in 



In the previous sections we have tried to define and describe equilibrium and driven interf 
general terms as briefly as possible. There are extensive review articles [13] and books 
on this subject and the reader is referred to them for details. In this thesis, we study a special 
class of interfaces which are produced by non-uniform fields which vary sharply in the direction 
perpendicular to the interface. The primary function of this field is to break the translational 
symmetry along the growth direction. Breaking this symmetry implies a non zero energy cost for 
long wavelength fluctuations of the interface which primarily tends to stabilize the sort of planar 
interfaces discussed in Section 1.3. Surprisingly, however, we discover, in each case, the existence 
of residual fluctuations which produce novel phenomena. 

In Chapter 2, we specialize our discussion of interfaces to interfaces in the Ising model. We 
introduce the model and discuss the equilibrium phases. The importance of Ising model in the 
sense that a variety of other statistical mechanical systems can be described by it, is then pointed 
out. In particular, reference is made to the mapping to lattice gas and binary alloy. Next we 
review existing literature on interfaces in the ferromagnetic Ising model. A study of the physics 
of single surfaces or interfaces begins with the characterization of the shape of the interface. 
In this regard, we study the structure of the Ising interface. We derive the interfacial profile 
using the continuum 0^ theory. We discuss the stability of a flat interface in the two-dimensional 
Ising model. The role of fluctuations is pointed out in detail. The time-dependent properties of 
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an inclined interface separating up and down spin regions in a two-dimensional nearest-neighbor 
Ising model evolving under Glauber dynamics in a non-zero but uniform field is discussed. We 
show that the for a constant magnetic field, the equation of motion for the interface reduces to 
the KPZ equation as expected. In the limit of large exchange coupling, the model reduces to 
the single-step model for ballistic growth and thence to the asymmetric exclusion process which 
describes a driven diffusive system of hard core particles on an one- dimensional lattice. The drift 
velocity of the interface is found as a function of field, temperature and inclination, and interface 
correlation functions are related to sliding tag correlation functions in the particle system. The 
non-linear dependence of the current on density leads to kinematic waves which involves moving 
density fluctuations. 

In Chapter 3, we study the steady state structure and dynamics of a two-dimensional Ising 
interface placed in an //7/7omogeneous external field translated with velocity Ve- The non-uniform 
field has a profile with a fixed shape which is designed to stabilize a flat interface. For small 
velocities the interface is stuck to the profile and is rippled with a periodicity which may be 
either commensurate or incommensurate with the lattice parameter of the square lattice. For a 
general orientation of the profile, the local slope of the interface locks in to one of infinitely many 
rational directions producing a "devil's staircase" structure. These "lock-in" or commensurate 
structures disappear as Ve increases through a kinetics driven commensurate - incommensurate 
transition. For large Ve the interface becomes detached from the field profile and coarsens with 
KPZ exponents. The complete phase -diagram and the multifractal spectrum corresponding to 
these structures are obtained numerically together with several analytic results concerning the 
dynamics of the rippled phases. The fact that interfacial fluctuations are partially suppressed 
by the external field is manifested in the exact agreement between a mean field theory and 
simulation results. However small interfacial fluctuations produce a dynamical phase diagram 
showing infinitely many dynamical phases and dynamic phase transitions. 

Next we turn our attention to more realistic systems viz., solid - liquid interfaces. In Chapter 
4, we introduce and describe atomic systems and discuss the need to study such systems. Before 
we go on to study interfaces in atomic systems we use this chapter to introduce the models and 
computational techniques (e.g. Monte Carlo and molecular dynamics simulations) that are used in 
this thesis. The relevance of these simulation techniques in the context of the thesis is pointed out. 
Specific techniques to simulate hard systems are also described. Bulk thermodynamic, structural 
and dynamic quantities that characterize different phases are described in some detail. We 
reproduce existing data in order to validate our computational methods to be used for subsequent 
calculations. 

In Chapter 5, we show how one can create a solid - liquid interface using a non - uniform 
external field. This is a direct analog of the Ising interface studied in Chapter 3. In contrast 
to the Ising interface, elastic deformations of the solid are allowed and play a prominent part in 
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determining tlie properties of the interface. In fact, we sliow tliat a tliin strip of solid placed within 
a bath of its own liquid relieves stress by novel interfacial fluctuations which involve addition or 
deletion of entire lattice layers of the crystal. Local analogues of the bulk quantities defined in 
Chapter 4 are calculated along the length of the system in the direction perpendicular to the 
interface. The "layering" transition is shown to be a generic feature of systems interacting via 
any interatomic potential. 

In Chapter 6, we show how these interfacial fluctuations influence mass, momentum and energy 
transport properties across the interface. Tiny momentum impulses are seen to produce shock 
waves which travel through the liquid solid interface and causes the spallation of crystal layers 
into the liquid. We show how kinetic and energetic constraints prevents the spallation of partial 
layers from the crystal, making this process suitable for creating nano-layer and coatings or for 
making nano-wires. We also study heat transport through the liquid solid interface and obtain 
the contact or Kapitza resistance of the interface as a function of the depth of the potential 
well. 



2 Interfaces in the Ising model 



After an introduction to the structure and properties of interfaces in general, in this chapter we 
specialize our discussion to a study of interfaces in the particularly simple context of the Ising 
model. The Ising model [izi tries to imitate behaviour in which individual elements (e.g., atoms, 
animals, protein folds, biological membrane, social behavior, etc.) modify their behavior so as to 
conform to the behavior of other individuals in their vicinity. The Ising model has more recently 
been used to model phase separation in binary alloys [17] and spin glasses jisf. In biology, it can 
model neural networks flocking birds 0|, or beating heart cells 0|. It can also be applied 
in sociology 

In this chapter we start by giving an introduction of the Ising model and its importance in 
explaining several phenomena. We discuss briefly the mapping of the Ising model to binary alloy 
and lattice gas models. We then discuss the structure and stability of the Ising interface in the 
low temperature limit when the interface may be described by a single valued "height" function 
h{x, t). Finally we discuss in detail the effect of a homogeneous external field on the Ising interface 
with special emphasis on the nature of fluctuations of the interface about the mean position h{t) 
of the interface. 



2.1 Ising model 

The Ising model was proposed in the doctoral thesis of Ernst Ising, a student of W. Lenz. Ising 
tried to explain certain empirically observed facts about ferromagnetic materials using a model 
proposed by Lenz (2^. It was referred to in Heisenberg's j3| paper which used the exchange 
mechanism to describe ferromagnetism. The name became well-established with the publication 
of a paper by Peierls [26], which gave a non-rigorous proof that spontaneous magnetization must 
exist. A breakthrough occurred when it was shown that a matrix formulation of the model allows 
the partition function to be related to the largest eigenvalue of the matrix (Kramers and Wannier 
IztI ]. Montroll 0: Kubo js^]. Kramers and Wannier calculated the Curie temperature 
using a two-dimensional Ising model, and a complete analytic solution was subsequently given by 
Onsager j3l| . 

Consider a lattice in d dimensions of sites i labelled 1,2, A/'(i7), which we will take to 

be hypercubic, unless otherwise stated. The degrees of freedom are classical spin variables. Si, 
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residing on the vertices of the lattice, which take only two values : up or down, or more usefully. 

Si = ±1 (2.1) 

The total number of states of the system is 2-^'^^^ The spins interact with an external field 
(in principle varying from site-to-site) Hi and with each other through exchange interactions 

Jij,Kijk, .. which couple two spins, three spins, etc. 

A general form of the Hamiltonian is 

— Ti-a = ^ HiSi + ^ JijSiSj + ^ KijkSiSjSk + ... (2.2) 

ien ij ijk 

The free energy is given by 

J^n(r, Hi, Jij...) = -keTlogTre-^^^ (2.3) 

where the trace operation is a sum over all possible configurations of the spin Si, and thermody- 
namic properties can be obtained by differentiating the free energy. 

To discuss the phase transition of the Ising model, we consider the simple case of the nearest 
neighbour Hamiltonian 

-nn^Hj2Si + jJ2 Si^J (2-4) 

i=l <ij> 

where we have assumed that the external magnetic field H is uniform in space, and that the only 
interaction between spins is that between neighbouring spins, with strength denoted by J. We 
consider the case J > 0. With a uniform external magnetic field, we can define the magnetization 
or magnetic moment per site, M : 

Now consider the case when H — 0. The spins then seek configurations that at constant T, 
will minimize the Helmholtz free energy J^q — £ — TS. At high T, J^ci is clearly minimized by 
maximizing the entropy, S. The maximally disordered state has the highest entropy, implying that 
the equilibrium state at high temperatures is the paramagnetic state with no average alignment 
of spins, i.e., no magnetization. At low temperatures, the internal energy dominates over TS, 
and the state that minimizes J-'n is one that minimizes £. States which minimize S have a 
non-vanishing magnetization, and the low temperature equilibrium phase is the ferromagnetic 
phase with nonzero magnetization M. At some temperature %, there is a phase transition 
from the entropy dominated paramagnetic state to the energy dominated ferromagnetic state. 
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The magnetization M is called the order parameter of the ferromagnetic phase and takes two 
degenerate values ±M(T) ^ ±1 as T ^ 0. The two degenerate states correspond to cases 
where the spins are either all up {Si = +1) or all down [Si = —1). At T = 0, all spins are 
perfectly aligned to either of the two possible values. At any nonzero temperature, domains of 
one phase in the other are stabilized due to entropy. These domains cost positive interfacial 
energy and grow in response to external fields. 

2.1.1 Mapping to other models 

The usefulness of the Ising model stems to large extent, from its generality. The Ising model 
maps onto a large number of simple statistical models for complex phenomena in material science 
as well as biology, economics and sociology. We show below three common mappings which are 
included as illustrative examples jljj. 




A. Ising antiferromagnetic model : If we replace J in Eqn. 2.4 by —J (J > 0) then energy is 
lowered for neighbouring spins with different signs. At low temperatures neighbouring spins are 
antiparallel. This becomes the antiferromagnetic model. In bipartite lattices where the lattice 
naturally breaks up into two interpenetrating sub lattices (e.g., square lattice in two dimensions 
and the body centered cubic lattice in three), the Ising antiferromagnet is mathematically identical 
to the ferromagnetic case. Indeed for, = 0, if we write every alternate Si as —Si, then Hn is 
the same as before because the extra minus sign cancels that from — J. In other words, if if = 0, 
the sign of J does not affect the thermodynamic potential. The order parameter for the Ising 
antiferromagnet is the "staggered" magnetization, Ms, which is the difference of the individual 
magnetizations of the two sub lattices. The analogue of the external field, H, is the "staggered 
field", Hg, which has opposite signs for each sublattice. 

B. Model of binary alloy : Let Si = 1 represent an atom of type A at point i and Si = —1 
represent an atom of type B at point i. Let gaa, ^bb, ^ab be the interaction energies between 
neighbouring atoms. A binary substitutional alloy which has either A or B types of atoms at 
each site may be described by the Ising model [3z]. The energy of a pair of neighbouring atoms 
i,j can be written as 




JSiSj - ):H{S, + Sj) + /C 




(2.6) 



where C, is the number of nearest neighbours for each spin. 
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Comparing the energies of a two atom cluster containing either AA, BB or AB types and 
solving for J, H and /C we get : 

1 1, , 

J — -<^AB — -{<^AA + '^BB) 

A- = -Cab + ^(eAA + esBj 

H = ^{eBB-eAA) (2.7) 

Using Eqn. 2.7, Eqn. 2.4 now describes the binary alloy. If {eAA + eBB)/2 < tAB{J > 0), at 
sufficiently low temperatures the atoms A will gather together while atoms B will form another 
cisuter. Hence we have a model for a de-mixing transition. If J < and the temperature is 
low, the atoms A and B will arrange alternately, leading to a order-disorder transition. The 
paramagnetic phase at high temperature of course, represents a randomly substituted alloy - the 
mixed phase. 



C. Lattice gas model : Let Si = 1 represent a molecule at site i, and Si = —1 represent no 
molecule at site i. We assume that the same site cannot accommodate two or more molecules. 
Let the interaction energy of neighbouring molecules be — e. Then we may write the energy 
between two neighbouring lattice sites as 

- J - — + /C = -e , S, = Sj = l 



c 

-J+^ + /C = , Si = Sj = l 

J + K. = , Si = l, Sj = -1 (2.8) 



Solving these, we get 



c 

— ( 

4 



H=^e (2.9) 



Therefore, Eqn. 2.4 is also a model of gas molecules, where — e represents the attraction between 
the gas molecules. The restriction of at most one molecule per site represents short-range 
repulsion. The high density state, i.e. {Si) > corresponds to the liquid phase while the low 
density state, i.e. (Si) > corresponds to the gas phase. Although this model may seem rather 
contrived since there is no underlying lattice in a real liquid-gas transition, surprisingly, the critical 
properties of this model turn out to be identical to real experimental gases due to universality 

0,3- 
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2.1.2 Structure of interfaces in the Ising model 

As mentioned above, at any non-zero temperature the Ising model contains domains populated 
by either one of the two ground states. At low enough temperatures, these domains are separated 
by well defined interfaces which fluctuate and move as one phase grows at the expense of the 
other. In this section, we shall study Ising interfaces for zero external field, where the two phases 
on either side of the interface are degenerate in energy. We shall also implicitly assume here and 
elsewhere in this thesis that the temperature is low enough compared to the coupling parameter 
J so that the interface can be represented either as a flat boundary or at least as a single valued 
curve which does not loop back on itself (i.e., contains no overhangs). To find the structure and 
behavior of the Ising interface we start with the Ginzburg-Landau free energy expansion 

T = j dv[-^<P{vf + ^0(r)^ + ||V0(r)p] (2.10) 

where the coefficients \l/, c and B are related to the microscopic parameters and is the space 
dependent coarse grained magnetisation. In particular ^ vanishes at the critical point given by 
%. The function (/)(r) which minimizes JF is determined by the equation : 

J^.?/_A^ = (211) 

where (p^ = dcp/dvi. The repeated index i and = {x,y,z) are summed over in Eqn. 2.11. 
The resulting equation is 

-^0 + c0^-5V^0 = O (2.12) 

with the boundary conditions that the system is uniform away from the interface. Far away from 
the interface, we expect the two phases to have their equilibrium values of the composition. 
Assuming that is small we get, 

00 = ±^ (2.13) 

Consider now a one-dimensional concentration variation, (p{z), with boundary conditions d(f)/dz - 
at 2; = ±00, which ensure the equilibrium phases at these limits and represents a flat interface. 
The solution is then 

\& z 

-tanh- (2.14) 



where the width of the interfacial region is given by ,^ = y^2Bj^, which is the bulk correlation 
length and diverges as \& — > or T ^ 7^. Thus, for z — > ±00, approaches its equilibrium, 
uniform values of ±^/¥/^ as shown in the Fig. 12.11 The width of the interfacial region is 
proportional to 1/v^, so that as the critical point is approached (\E' — > 0) the interfacial width 
diverges. 
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Figure 2.1: A plot of the interfacial profile 
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(b) 



Figure 2.2: (a) Boundary line for a square lattice, (b) An enlargement of a portion of (a). 



2.1.3 Stability of the Ising interface 

To discuss the stability of the flat Ising interface described above, let us consider a square lattice 
of Ising spins with a boundary line (interface) separating up and down spin regions as shown in 
the Fig. 12.21 The energy of the interface is its total length times 2 J. The total length is 



L = L^ + Y,\yk\ (2.15) 

k=l 

Here, yt is the length of the k — th boundary counting from the left. In Fig. 12.21 -^i = 2,y2 = 
l,?/3 = 0,?/4 = — 2, ?/5 = —1,.... etc. Therefore the thermodynamic potential JF of the full 
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interface in the limit of large Lq, is 

JF = —ksTlnZ 

oo oo oo 

2 = J2 Yl E (2.16) 

i/i=-oo j/2=-oo yLQ=-oo 

Now, in the k-th horizontal position, the height of the boundary is 

k 

hk = Y,yj (2-17) 

i=i 

Hence, the difference in height between two points on the interface is 

k+n 

Ah = hk+n -hk= ^ Vj (2.18) 
j=k+i 

If n « Lq, each will be independent variables. Using the central limit theorem we get the 
distribution of A/i : 

\j2Txa 
= n{ip) 

n 

2 sinh'^ {J /keT) ^^'^^^ 

Hence the interface fluctuates violently in the same manner as a one dimensional random walker. 
The amplitude of fluctuation in the middle is about a/A)/ sinh( J/Ze^T). Although ^/L^ is much 
smaller than Lq, it is still not a microscopic length scale. Therefore, this interface is a wiggly 
and rough curve, not smooth and flat. In this analysis we have not considered "overhangs" or 
isolated regions. However the conclusion that the interface is rough is valid even if we include 
these effects. 



2.2 Driven Interfaces in the Ising model 

When there is an external field, the interface in an Ising model moves in a way that increases 
the size of the domains which are aligned in the direction of the field. The general problem of 
such a driven interface even for a simple model system like the Ising model is quite complex. In 
this section we show how, by using simplifying assumptions one can study such a driven interface 
in the two-dimensional Ising system. This is a reproduction of earlier work on constant field 
driven interface by Majumdar and Barma 0. E^. E^. E^. 'n the limit that the field H and 
temperature T are both much smaller than the nearest neighbour exchange coupling J (assumed 
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isotropic), one can neglect overhangs so that the interface can as before be represented as a 
single valued function. The interface moves in steady state with a velocity Vf which depends 
on the external field H and on the orientation 9 of the interface, measured with respect to the 
horizontal (Fig. 12. 3j) . As the Ising spins obey single-spin flip Glauber dynamics, and H,T « J 
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Figure 2.3: An Ising interface (bold line) between up and down Ising spins. Since J is very large, 
only spins at the corners, like the circled ones (coloured red) can flip. The orientation Q of the interface 
is also shown in the figure. 

holds, the interface moves by "corner flips" where a vertical bond and a horizontal bond at a 
corner exchange places. 

We shall show how such driven interfaces may be described using the KPZ equation. We will 
derive this dynamical equation first from the time dependent Ginzberg-Landau equation (TDGL) 
jsT], and then once again after mapping the corner flip dynamics of the two-dimensional Ising 
model to the dynamics of a one-dimensional gas of hard core particles (the exclusion process) in 
the limit H/J,T /J 0. We shall use the latter mapping to deduce some important dynamical 
properties of this interface. 

2.2.1 Derivation of KPZ from TDGL 

We consider, for the moment, that the inclination of the interface 9 = 0. The interface is given by 
the function h{x,t). The coarse grained, space dependent magnetization (p is given by a function 
which is uniform everywhere except near the interface, h(x,t) such that = (p^y — h{x,t)). The 
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function (p has a profile similar to Fig. 2.1. Model A dynamics [34] for </> then implies 



where 



j rfr[ai02 + a20^ + a3(V0)' - (2.21) 



is the coarse-grained Hamiltonian (which, for the moment, ignores the lattice) of an Ising system 
in a external field H and is a Gaussian white noise with zero mean and 

< C(r,t)C(r',t') > = 2kBTV5{T -Y')5{t-t') (2.22) 

Using Ht in Eqn. l2.20l taking y — h{x, t) = u and converting all derivatives to derivatives over 
the profile h{x,t), we have, 

dh d'^h 
-0'(^^)^= - r[2ai0(z/)+4a20=^(z/)-2a30"(z/)+2a30'(z/)^ 

dh 

- 2a30»(— )^-iJ]+C(r,t) (2.23) 

where primes denotes derivatives with respect to v. We then choose a (p dependent mobility 
r contributing to the lowest order in (p consistent with symmetry viz. F = Fq + Fi(V0)^. 
Substituting for F and integrating both sides of the equation with respect to h between limits 
(h — x/2) and (/i + x/2) i e. over the interfacial region, we finally get an equation of motion for 
the interface. 

dh , S^h , fdh\-i 



where Ai,A2 and A3 are parameters which involves the integral of 0(i/),0'(z/) and 0'(z/) over the 
regions of the interface. Note that Eqn. 12.241 lacks Galilean invariance-^ h' ^ h + ex, x' 
X — A2et, t' t. Eqn. 2.25 is the familiar KPZ equation and therefore the interface grows with 
KPZ exponents a = 1/2,13= 1/3. 



2.2.2 Mapping to the Exclusion Process 



In the exclusion process we consider a one-dimensional lattice of A/'s sites, p 
are occupied by particles, and assume periodic boundary conditions jssi, 
exclusion process, a particle chosen at random attempts to hop with probability p to the right 



which Afp = pMs 
3?! . In the simple 



^Note that f{Y,t) in Eq. H2.24|l generates a space and time dependent, (annealed) random, Galilean boost. 
Random Galilean transformations often lead to multifractal steady states; see for eg. U. Frish, Turbulence 
(Cambridge University Press, 1995) 
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and probability q to the left, with p + q = 1. Because of the hard core constraint, the hop actually 
takes place if the sought site is vacant.A/^ such attempted hops constitute a single time step. 
In the steady state, every configuration of of the Mp particles is equally likely ISl -Let 

us label the particles n = 1,2, ,Afp sequentially at t = 0. The ordering is preserved by 

the dynamics of the exclusion process. The configuration of the system is specified by the set 
{h(n)} where h(n) denotes the location of the nth particle. The corresponding interface model, 
called the particle height (PH) model, is obtained by interpreting the tag label n as a horizontal 
coordinate, and h(n) as a local height. Each configuration {h{n)} then defines a one- dimensional 
interface in the form of a staircase inclined to the horizontal with mean slope tan^ = 1/p- The 
interface coordinates satisfy h{n + l) > h{n)-\-l, and the periodic boundary conditions translate 
into h{n + J\fp) = h{n) iA/'^.The evolution rule is as follows: in each time step, h{n) tends to 
increase (or decrease) by 1 with probability p (or q); it actually increases (or decreases) if and 
only if h{n + 1) — h{n) > 1 (or h(n) — h{n — 1) > l).ln the unbiased case {p = q = 1/2). the 
interface does not move with a net velocity, but fluctuates around its initial position. But in 
the biased case {p ^ q), the interface moves vertically with the particle drift velocity Vp. The 
probabilities p and q can be related to the Ising model parameters on noting that ratio of the 
rates of up-down and down-up flips is exp(— 2/3if); thus p = exp(/5if )/(2 cosh(/?if) and the 
bias A = p — q = tanh(/9if). 

Vp =< {h{n, t) - h{n, 0) >= (1 - p)(p - q) (2.25) 

where h(n,t) is the location of the nth particle at time t. 



2.2.3 Some Exact Results 

It is straight-forward to derive 0, Q the expression for Vp and the particle current j. The 
particle velocity, 

Vp =< ih{n, t) - h{n, 0) >= (1 -p)ip- q) (2.26) 

which is the statement of the fact that in order for a particle to move one needs a bias and a hole 
to move into. Although we have defined the exclusion process in terms of particle hopping, it can 
equally well be viewed as the backward motion of holes(sites on which there are no particles). 
The drift velocity for labelled holes is Vho = —{p — q)p and the steady state current is 

j = = (1 - p)Vho = p(l - p)(p - q) (2.27) 

In time t, the mean vertical shift of the interface is given by the average number of particles 
which pass by a given hole, namely, p{vp — Vho)t- The magnitude of the normal velocity is then 

Vf = p{vp — Vha)cos6 = j{sm6 + cos6) which may be rewritten as 

tanh i3H ^ 
~ (sec 9 + cosec6') ' ^ ' ^ 
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Figure 2.4: Plot of j versus density 



The space dependent coarse-grained density p{n) and current J(n) are related through the 
continuity equation 

dp dJ{n) 



dt ~'~ dn 

Now the current J(n) can be written as |4l| . 



0. 



Jin) 



(2.29) 



(2.30) 



where the first term is the usual diffusion of particles with diffusion constant D, t] \s a Gaussian 
white noise and j(p) is the current due to the local density p = p{n). We may replace j(p) by its 
macroscopic value (Eqn. 2.27) and expand it in power series in the density deviation A = p — p 
around the average density p = Mp/M. 



Kp) = J^JmA" 



(2.31) 



with m\jm = d'^J/dp'^. The density fluctuation field A{n,t) is related to the height variable 
h{n,t) which measures the transverse displacement of the interface: 



h{n, t) 



A{n',t)dn'. 



(2.32) 



no 



Differentiating Eqn. 2.32 and using Eqn. 2.29 and Eqn. 2.30 one can show that h satisfies the 
equation by KPZ The equation of motion of the interface is 



dh 
dt 



ori'^ ^-^ on 

m 



+ V{n,t) 



(2.33) 
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Interfaces in the Ising model 



We consider the height fluctuation correlation function S{n,t) defined by 

S^{n,t) = < {h{n,t) - h{n,0) -VptY > (2.34) 

The large-distance large-time behaviour of S(n,t) depends strongly on the first few coefficients 
in Eqn. 2.33. The constant jo is equal to the macroscopic current (Eqn. 2.27) and determines 
the mean rate of growth of the interface. It can be eliminated by the boost h h + jot. 
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Figure 2.5: The second moment of the height distribution S'^/t as a function of time t. 

The first-order gradient term ji is equal to U, the speed of local density variations or kinematic 
waves j4^. Kinematic waves, which transport density fluctuations in the lattice system, corre- 
spond to the movement of transverse fluctuations in the interface problem and owes its existence 
to the conservation of particle number. The term ji can be eliminated from Eqn. 2.34 by making 
the Galilean shift x — > x + jit. This has the effect of moving to a coordinate system in which the 
wave is stationary. If the jo term is eliminated but ji is not, it can be shown that S{0,t) grows 
as t^/^. Since the sliding of the wave is the dominant physical effect so S{0,t) is approximated 
by the equal time correlation function S{n = Ut,0). Since equal-time height fluctuations are 
determined by fluctuations of particle number, these grow as n^^^, and it follows that 

S{0,t) ^ {Uty/\ (2.35) 

We have performed Monte Carlo simulations with random sequential updates of the 

particles to reproduce these exact results. This constituted the benchmarking and validation for 
our codes which were subsequently used to study the system discussed in Chapter 3. In Fig. 
2.4 we show a plot of the measured current as a function of the density together with the exact 
result. Ms = 10000 was used in this run. In Fig. 2.5 we have shown our results for the second 
moment of the height distribution function. Our Monte Carlo results clearly depicted the t^/^ 
behaviour. We used J\fs = 80000 for this calculation. 



3 Ising Interfaces in spatially varying 
external field : Novel fluctuations 



In the last chapter we introduced the Ising model and the Ising interface. We studied the structure 
of this interface in two-dimensions and discovered that at any non zero temperature the interface 
is never flat but fluctuates rather strongly. The random fluctuations of the interface diverge as 
the system size as L^/^ and with time as t^/^. 

In this chapter, we discover how to control these interfacial fluctuations and create an essentially 
flat interface. The ability to grow flat solid surfaces is often of major technological 

concern, for example, in the fabrication of magnetic materials for recording devices where surface 
roughness jl^l causes a sharp deterioration of magnetic properties. We show that non-uniform 
external fields-^ which have a sharp gradient at the position of the interface, suppresses height 
fluctuations to a large extent. This is a direct consequence of the fact that such spatially 
varying fields break the translational symmetry of the interface in the direction of growth so 
that long wavelength fluctuations now do cost non zero interfacial energy. Surprisingly however 
we show that this suppression of long wavelength fluctuations generate new short wavelength 
structures and patterns. Driving the ineterface by using a moving field profile causes interesting 
dynamical transitions between patterns which are either commensurate or incommensurate with 
the underlying lattice structure over which the Ising model is defined. 

Commensurate-in -commensurate (C-l) transitions have been extensively studied over al- 
most half a century following early experiments on noble gases adsorbed on a crystalline substrate 
I49I 1 eg. Kr on graphite. Depending on coverage and temperature, adsorbates may show high 
density periodic structures the reciprocal lattice vectors (RLVs) of which are either a rational 
(commensurate) or irrational (in -commensurate) multiple of a substrate RLV. By changing ex- 
ternal parameters (eg. temperature) one may induce phase-transitions between these structures. 
Recently, the upsurge of interest in the fabrication of nano-devices have meant a renewed interest 
in this field following a large number of experimental observations on "self-assembled" domain 
patterns (stripes or droplets) on epitaxially grown thin films for eg. Ag films on Ru(OOOl) or 

^ There are several practical examples where non-uniform fields drive interfaces. Some of them include 
zone purification of Si where the controlled motion of a temperature field profile is used to preferentially 
segregate impurities ,4S|, magnetization of a bar of iron with a permanent magnet, phase transitions 
induced by a travelling heat (welding) or pressure (metamorphosis of rocks) fronts etc. 
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Cu-Pb films on Cu(lll) [SC^ etc. The whole area of surface structure modifications and surface 
patterns has tremendous technological implications for example in opto-electronics, recording 
industry, coatings and paints, etc. 

Almost universally, C-l transitions may be understood using some version of the simple Frenkel - 
Kontorova 0| model, which models them as arising from a competition between the elastic en- 
ergy associated with the distortion of the adsorbate lattice and substrate -adsorbate interactions. 
A complicated phase diagram involving an infinity of phases corresponding to various possible 
commensuration ratios (rational fractions) is obtained as a function of the two energy scales. In- 
between two commensurate structures one obtains regions where the periodicity of the adsorbate 
lattice is in -commensurate. 

The structures obtained in our interfacial problem, on the other hand, undergo C-l transitions 
which are entirely dynamical in origin and can be controlled by adjusting the velocity of the 
interface. Short periodic structures are stabilized at lower velocities. As the velocity of the 
interface approaches a limiting velocity Voo, the patterns disappear and the interface begins to 
fluctuate over all length scales. Beyond Voo, KPZ behaviour takes over. 



3.1 The Model 

We show in Fig. 13.11 a one-dimensional interface h(x,t) between phases with magnetization, 
(f){x,y,t) > and (l){x,y,t) < 0, in a two-dimensional Ising model on a square lattice^ obeying 
single-spin flip Glauber dynamics in the limit H/J,T/J 0. An external non-uniform field 
is applied such that H = Hmax in the +ve and —Hmax in the -ve regions separated by a 
sharp edge. The edge of the field (i.e. where the field changes sign) lies at Se- The front or 
interface, h(x,t), (no overhangs I) separates up and down spin phases. The interface is shown 
as a bold curved line with the average position Sj. To move the interface we move the edge with 
velocity Ve', in response the front moves with velocity Vf. Parts of the front which leads (lags) the 
edge of the field experience a backward (forward) force pulling it towards the edge. The driving 
force therefore varies in both space and time and depends on the relative position of the front 
compared to that of the edge of the dragging field. In the low temperature limit the interface 
moves solely by random corner flips [10] , the fluctuations necessary for nucleating islands of the 
minority phase in any region being absent. We study the behaviour of the front velocity Vf and 
the structure of the interface as a function of and orientation. 

Naively, one would expect fluctuations of the interfacial coordinate h{x,t) to be completely 
suppressed in the presence of a field profile. This expectation, is only partially true. While, as we 
show below, a mean field theory gives the exact behaviour of the front velocity vj as a function 

'^(j) represents any scalar order parameter with h{y,t) a field conjugate to 0, eg. number density and 
chemical potential for a growing gas -solid surface. 
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Figure 3.1: An Ising interface h{x,t) (bold curved line) between regions of positive (marked +) 
and negative (marked — ) magnetization in an external, inhomogeneous field with a profile which is as 
shown(dashed line). The positions of the edge of the field profile and that of the interface are labelled 
Se and Sf respectively. 
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Figure 3.2: The interface velocity vj as a function of the velocity of the dragging edge for 
Ms = lOO(n), 1000(0), 10000(+) and p = 0.5. All the data (□,<>,+) collapse on the mean field 
solution (dashed line). Inset shows the graphical solution (circled) of the self-consistency equation for 
Vf] dashed line represents Vf = Vf. 



of Ve', snnall interfacial fluctuations produce a dynamical phase diagram showing infinitely many 
dynamical phases. For Ve < the interface is stuck to the profile Vf = v^- The stuck phase 
has a rich structure showing microscopic, "lock-in", commensurate ripples. These disappear at 
high velocities through a dynamical C-l transition. For large Ve the interface detaches from the 
profile and moves at velocity corresponding to an uniform field of magnitude H — H^ax 
independent of and coarsens with KPZ exponents. 
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3.2 Continuum Description 

In the continuum description of this model, we need to substitute H = H(h,t) = Hmaxf{h,t) 
with f{h,t) = ta,nh.{(h — Vet) / x) where x is the width of the profile (see Fig lS.!]) in the equation 
for the Hamiltonian in Eqn. 2.22. Carrying out an analogous procedure for obtaining the equation 
of motion for the interface we obtain a generalized KPZ equation for the dynamics of the interface 
in a spatially varying field. The result analogous to Eqn. 12.241 is 

^ = Ai^-A2(^) f{Kt)-X,f{h^t) + CM (3.1) 

A mean field calculation amounts to taking h = h(t) i.e. neglecting spatial fluctuations of the 
interface and noise. For large times (t oo), h — > Vft, where vj is obtained by solving the 
self-consistency equation; 

. fiVf — Ve)t\ 

Vf = lim — Aatann — 

■> t^oc \ X ^ 

= -Xssigni^Vf - Ve) (3.2) 

For small Ve the only solution to Eqn. 13.21 is Vf = and for Ve > foo, where v^o = A3 we get 
Vf = \3 = Voo- We thus have a sharp transition (Fig. I3.2jl from a region where the interface is 
stuck to the edge to one where it moves with a constant velocity. How is this result altered by 
including spatial fluctuations of /i ? We answer this question by mapping the interface problem 
to an asymmetric exclusion process 0, 0, '^^J and studying the dynamics both analytically and 
numerically using computer simulations. 



3.3 Beyond Mean Field Theory 

The mapping to the exclusion process follows [36, E3 0, 0, Q as in the previous chapter 
by distributing Mp particles among J\fs sites of a 1-d lattice. The particles are labelled i = 

1,2, ,J\fp sequentially at t = 0. Any configuration of the system is specified by the set of 

integers {hi} where hi denotes the location of the ith particle. In the interface picture i maps 
onto a horizontal coordinate {x in Fig. 13. Ij) . and hi as the local height h(x). Each configuration 
{hi} defines a one-dimensional interface inclined to the horizontal with mean slope tanOj = 1/p 
where p = Mp/Mg- The hi satisfy the hard core constraint /ij+i > hi + 1. The local slope 
near particle i is given by /ij+i — hi and is equal to the inverse local density pi measured in a 
region around the i^^ particle. Alternatively, one associates a vertical bond with a particle and a 
horizontal bond with a hole 00, in which case, again, we obtain an interface with a slope 
tan 6'j: = p/{l — p). The two mappings are distinct but equivalent. Periodic boundary conditions 
amount to setting /ij+at = hi ^Mg- Motion of the interface, by corner flips corresponds to the 
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hopping of particles. In each time step {J\fp attempted hops with particles chosen randomly and 
sequentially jso], ll^), hi tends to increase (or decrease) by 1 with probability p (or q); it actually 
increases (or decreases) if and only if hi+i — hi > 1. The dynamics involving random sequential 
updates is known to introduce the least amount of correlations among hi which enables one 
to derive exact analytic expressions for dynamical quantities using simple mean field arguments 
(see chapter 2). The right and left jump probabilities p and q {p + q = 1) themselves 
depend on the relative position of the interface hi and the edge of the field profile i/p + Vet. 
Note that this relative position is defined in a moving reference frame with velocity the 
instantaneous average particle velocity defined as the total number of particles moving right per 
time step. We use a bias Aj(t) = p — q = Asign(/;,j — i/p — vj) with A = 1 unless otherwise 
stated. Our model is thus a generalization of the one-dimensional exclusion process with space 
dependent, dynamic jump probabilities 0, Ezj- addition to the front velocity vj, we also 
examine the behaviour of the average position, 



< h{t) >-- 



hi{t) 



(3.3) 



and the width of the interface: 



= E < (^^W- < >)' > (3-4) 

i=l,Afp 

as a function of time and system size A/'s. Here, < hi(t) >= i/p + Vet. Angular brackets denotes 
an average over the realizations of the random noise. Note that the usual particle hole symmetry 
for an exclusion process jlo|, 0, 10 is violated since exchanging particles and holes changes the 
relative position of the interface compared to the edge. This violation is, of course, completely 
equivalent to the breaking of translation symmetry of the interface in the growth direction as 
mentioned before (chapter 1). 

We perform numerical simulations of the above model for A/'s upto 10^ to obtain vj for the 
steady state interface as a function of fg as shown in Fig. 13.21 A sharp dynamical transition from 
an initially stuck interface with vj = f e to a free, detached interface with vj = v^c = A(l — p) 
is clearly evident as predicted by mean field theory. The detached interface coarsens with KPZ 
exponents js^j. Note that, even though the mean field solution for Vf{ve) neglects the fluctuations 
present in our simulation, it is exact. The detailed nature of the stuck phase [vf = and w 
bounded) is, on the other hand, considerably more complicated than the mean field assumption 
h{x,t) = h{t). Below we analyse the nature of the stuck phase starting from the ground state 
configurations at v^, = 0. 
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3.3.1 The Ground State and the Devil's Staircase 

The ground state of the interface in the presence of a stationary (t>e = 0) field profile is obtained 
by minimizing E = YliA^i ~ "^1 P ^ c)^ with respect to the set {hi} and c. This maybe shown 
from Eqn 13.11 by neglecting terms containing spatial derivatives of h; the resulting equation of 
motion, for small deviations of h from the edge may be derived from the effective Hamiltonian 
E. The form of E leads to an infinite range, non-local, repulsive, interaction between particles 
in addition to hard core repulsion and the minimization is subject to the constraint that hi be an 
integer. For our system, the result for the energy may be obtained exactly for density p = m/n, 
an arbitrary rational fraction. For even m we have the following expression. 



E 



1 

m 



--cr+{-+cY+{--cr+{-+cr 

m m m m 



+ ....+ (A_c)2 + (A + c)2 + (^-c)2 + c2 
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(3.5) 

m m 4 

Minimising with respect to c, we get c = Substituting for c we get the expression for 
minimized energy as 



6 2m m Am 4m^ 
Similarly for odd m we have the following expression for energy. 



(3.6) 
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(3.7) 



which when minimised with respect to c gives c = 0. Hence, the minimised energy for odd m is 
given by 



E = 1^(1 -A) (3-8) 
12 m^' 

The resulting ground state profiles are shown in Fig. 13.41 The lower bound for E(p) is zero which 
is the energy for all p = l/n. For irrational p the energy is given by \ira,m^oo E (m/n) = 1/12 for 
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Figure 3.3: Plot of E{p) for a Ms = 420 site system. 

both even and odd m which constitutes an upper bound (Fig. IS.Sjl . For an arbitrary < p < 1 
the system {{hi}) therefore prefers to distort, conforming within local regions, to the nearest 
low-lying rational slope 1/p interspersed with "discommensurations" of density = \p — p\ and 
sign +ve (-ve) if these regions are shifted towards (away) from each other by 1. A plot of p{p) 
shows a "Devil's staircase" structure [55, 56J. We observe this in our simulations by analysing the 
instantaneous distribution of the local density of the particle-hole system to obtain weights for 
various rational fractions. A time average of these weights then give us the most probable density 
p - distinct from the average p which is constrained to be the inverse slope of the interface. 

3.3.2 Dynamics of interfacial patterns for small Ve 

What are the low lying excitations of the ground state structures ? For low velocities and density 
where correlation effects due to the hard core constraint are negligible, the dynamics of the 
interface may also be obtained exactly. Under these circumstances the Mp particle probability 
distribution for the hi's, P{hi, h2, - ■ ■ , hj^f^) factorizes into single particle terms P{hi). Knowing 
the time development of P{hi) and the ground state structure the motion of the interface at 
subsequent times may be trivially computed as a sum of single particle motions. A single particle 
(with say index i) moves with the bias Aj(fet) which, in general, may change sign at /i < 
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Figure 3.4: Variation of < h{t) > with t for Ve = 0.025 and p = 1.0. Lines denote analytic results 
while points denote Monte Carlo data for p = 1/5 (uppermost curve), 2/5 and an incommensurate p 
near 1/3. Inset (a)-(c) shows the corresponding ground state interfaces {hi —i/p). The arrows in (c) 
mark the positions of two discommensurations. 

i/ p + Vet < h + 1. Then P{hi) satisfies the following set of master equations, 

P{hi) = -P{hi) + P{hi + 1) forhi>h + l 

P{hi) = P{K - 1) - + P{hi + 1) for hi = h, h + 1 

p[hi) = -P{hi) + P{hi - 1) for hi < h. (3.9) 

Note that the average position of the particle is given simply by < hi{t) >= Yl'h',=~oo hiP{hi) 
and the spread by w'^{t) = Yl'h =-00^^^" ^ ^i(^) >Y P{hi)- Solving the appropriate set of 
master equations we obtain, for « 1 the rather obvious steady state solution P{hi) = 
l/2{Shi,h+Shi,h+i) and the particle oscillates between h and h+1. Subsequently, when i/p+Ve t > 
h+1, the particle jumps to the next position and P{hi) relaxes exponentially with a time constant 
r = 1 to it's new value with h ^ h + \. For p = 1/n the entire interface moves as a single 
particle and the average position advances in steps with a periodicity of l/v^ (see Fig. I3.4|l 

In general, for rational p = m/n, the motion of the interface is composed of the independent 
motions of m particles each separated by a time lag of tl = 1/mVe- The result of the analytic 
calculation for small Ve and p has been compared to those from simulations in Fig. 13.41 for 
p = 1/5 and 2/5. For a general irrational p < 1/2, m ^ 00 consequently, 0. The hi's 

are distributed uniformly around the mean implying w"^ = 1/3 independent of system size and 
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time. For p > 1/2 the width w"^ = (1 — p)/3p since the number of mobile particles decreases by 
a factor of (1 — p)/p. 

In chapter 2, we encountered kinematic waves which were travelling density modulations in the 
exclusion process which, in the interface picture, represent height fluctuations which travel along 
the interface with velocity U = dJ/dp. In this case too, the forward motion of the interface is 
accompanied by the motion of discommensurations along the interface. The velocity of these 
kinematic waves in this system, however, is constrained to be equal to v^. 

As the velocity t>e is increased, the time lag, r^, decreases and the patterns begin to get 
distorted becomes comparable with the MC time step. The instantaneous values of p begins 
to make excursions to other nearby fractions and eventually becomes free. Since tl is smaller 
for large m, higher order fractions becomes unstable earlier. In Fig. 13.51 thus effect is clearly 
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Figure 3.5: Devil's staircase structure for two different velocities. Steps disappear as velocity of 
the edge increases. 

visible; the Devil's staircase for the p{p) curve for the higher velocity has fewer steps. Steps 
corresponding to p = m/n disappear (i.e. p ^ p) sequentially in order of decreasing m and the 
interface losses the ripple patterns. 

3.3.3 Dynamical transitions and the dynamical phase diagram 

The beginning and end points of each step in the curve (Fig. 13.511 mark the stability limits 
for each of the patterns which may be regarded as distinct dynamical phases. The locus of these 
points as a function of v^. trace out the dynamical phase diagram which is shown in Fig. 13.61 
As expected, fractions with higher values of m disappear as increases making way for simpler 
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Figure 3.6: The dynamical phase diagram in Ve and p plane. The numbers on the p axis mark the 
fractions p, which determines the orientation of the lock-in phase. The three regions white, black and 
grey correspond to the rippled, the disordered and the detached phases respectively. 

structures. We thus have true dynamical transitions from more complicated to simpler patterns 
as Ve approaches foo- Beyond Voo, the interface detaches and we get back the KPZ interface. 

The sequential smoothening of the Devil's staircase may be naturally described using the 
language of multifractals. For every Ve, we first obtain the generalized dimension Dk- For this 
purpose we obtain the scales /j which are the differences in p between subsequent rational fractions 
in the Farey construction (e.g. p = 0,1,1/2,1/3,2/3, etc). For every scale, the corresponding 
measure tt^ is given by the jump in p which depends on Ve and the structure of the Devil's 
staircase. The generalized dimension Dk is obtained for any given value of k by solving 

J2(^^/it'^''n = 1- (3.10) 
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Once Dk is obtained as a function of k, f{a), the multifractal spectrum of singularities, is 
obtained by a Legendre transform of {k — l)Dk (see Fig. 13. 7p . 

a{k) = ^[{k - 1)D,] 
f{a) = k-^[{k - l)Dk] -{k- 1)D, (3.11) 

In the vicinity of the primary steps corresponding to p = 1/n, p ~ {p — PmaxY where p^ax is the 
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Figure 3.7: The multifractal spectrum /(a) for two different velocities We/(1 — p) = 0.1 and 0.5. . 
Note that with increasing velocity, amin (see text). 

largest density at which the step is stable. This universal critical exponent ^ = D^o has been 
determined to be 0.71 ± .001 from our data at fe/(l — P) = 0.1. The exponent ^ determines how 
strongly p diverges away from its value at a step for small changes in p. Since the density p is 
related to the average orientation of the interface which is in turn determined by the orientation 
of the external field profile, ^ determines the stability of the ripple pattern to small changes in 
the external field. As Ve increases, ^ — implying that the 1/n steps become extremely stable 
at the expense of those corresponding to higher order rational fractions. 

The nature of the C-l transitions discussed here is essentially new and fundamentally different 
from that described within the Frenkel-Kontorowa (F-K) model. Indeed the control parameter 
for our phase diagram is t>e - a dynamical parameter. Also, the non local E{p) and the non linear 
constraint hi = integer makes it impossible to devise a natural mapping of this problem to an 
effective F-K model. We therefore use a different approach as discussed below. 
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Figure 3.8: A surface plot of p vs p and \fTK showing steps for the fractions 1/4,2/7 and 1/3. 
Note that the step corresponding to 2/7 vanishes at \/Tk > 0.5. 

3.4 Langevin Dynamics 

While the stability of mismatch domains [56J is decided, mainly, by competition between mechan- 
ical, long-ranged (elastic) and short-ranged (atomistic) interactions [55], dynamical ripples vanish 
with increasing t>e through increased fluctuations. We argue that it is sufficient to project the 
entire configuration space hi of the stuck interface onto the single variable p. As is obvious from 
the energy versus rho diagram (Fig. 13.311 -E'(p) has a structure similar to the free energy surface 
of a 1-d "trap" model [SSj often used to describe glassy dynamics. The distinction, of course, 
is the fact that the energies of the traps in this case are highly correlated. We now attempt to 
describe the dynamics of the stuck interface as the Langevin dynamics viz. 

, dF 

dp' 

<Vp'{i)Vp'ii') > = 2T5{t-t') (3.12) 
of a single particle with coordinate p' diffusing on a energy surface given by, 

Fip')=Eip') + Kip'-pr. (3.13) 



The particle is kicked by a Gaussian white noise {rip) of strength T. The second term, containing 
the modulus n, ensures that for infinitely large times p' always relaxes to the true global minimum 
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p' = p. At intermediate times, however, the system may get trapped indefinitely in some nearby 
low-lying minimum with p' = p if the noise strength T is not large enough. When the time taken 
for jumping between minima exceeds the residence time in any minimum because of enhanced 
noise, we have a fluctuation induced C-l transition (Fig. 13. 8j) from that density p. To show this 
we obtain the limiting value of p' averaged over realisations of the noise as a function of p and 
T from a numerical solution of Eqn. 13.121 While it is difficult to make a direct quantitative 
connection between T and v^, symmetry considerations would require T oc v^. In Fig. 13.81 we 
have plotted p{p,T) for 1/4 < p < 1/3 showing prominent steps for p corresponding to the 
rational fractions 1/4,1/3 and 2/7. The steps for 2/7 vanish at the square root of the noise 
amplitudes as 1/2 which is as expected. 

3.5 Behavior at the transition point 

We want to determine scaling form for w{t) at the transition point viz. the growth exponent P, 
the roughness exponent a and the dynamic exponent z. In the detached phase we know from 
renormalization group analysis that the exponents are in the KPZ universality class 0, viz. 
P = 1/3, a = 1/2 and z = a/P = 3/2. To determine these exponents at the transition point 
we make use of Family-Vicsek scaling relation |lo| 



Fig. 13.91 shows the variation of t/M^ with w{L,t)/M^ for different p and different system sizes 
Ms. The curves collapse onto one curve once an intrinsic width Wi, arising from finite-size and 
crossover effects [lOj, is subtracted out. The exponents were found to be KPZ. 

To understand why this happens we go back to our modified KPZ equation Eqn. 3.1 and make 
the transformation h = h' + vjt. We get 



Now substituting the mean field result for vj (Eqn. 13.211 . making use of the fact that at the 
transition point Vf = A3 and simplifying one can show that the above equation reduces to the 
familiar KPZ equation in h'. 



(3.14) 




(3.15) 
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Figure 3.9: Monte Carlo data {w - Wi)/M^ vs t/Ml for p = 1.0, p = 0.7 and with Ms = 1000 
[x,0] . 800 [0,o] , 400 [+,*] , 200 [0,A]. All the curves collapse to a single universal function 
showing KPZ scaling. 



po.3 



1 1 








.->■ 




1 1 





0.1 0.2 0.3 0.4 0.5 
P 



Figure 3.10: Devil's staircase structure for Ve = 0.05 and 1/x = 0.01, 1, 5. 

3.6 Overview, caveats and conclusion 

In this chapter we have shown how Ising interfaces in spatially varying, dynamical fields can give 
rise to interesting and new phenomena like dynamical C-l transitions. We hope that our study 
may lead to experimental work on real interfaces in carefully controlled and shaped external fields 
which may be used to carve static or dynamical patterns with technological applications. There 
are several problems which need to be sorted out before these techniques begin to have any 
practical application. 

Firstly, the sharpness of the field profile directly influences the "resolution" of the patterns, 
Even at small velocities, as the width of the profile x increases, the steps in the Devil's staircase 
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corresponding to particular patterns disappear (Fig. I3.10|) - simpler patterns being more stable as 
expected. In order to obtain the more interesting patterns therefore one needs the field profile to 
be atomistically sharp. Secondly, the most useful interfaces are two dimensional and therefore our 
studies need to be extended to higher dimensions and for other underlying crystal structures. We 
believe that there are exciting possibilities to be explored in this direction. Lastly, Ising interfaces 
miss two very important characteristics of real interfaces. (1) Elastic distortions and the effects 
of interfacial stress and (2) the possibility of particle transfers between the two phases. It is this 
last possibility that we explore in the subsequent chapters where we look at a two dimensional 
liquid solid interface stabilized by a spatially varying chemical potential field. 
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4 Atomic Systems 



4.1 Introduction 

In the last two chapters we have concentrated on interfaces in the two dimensional Ising model. 
Ising interfaces are particularly simple because they represent spatial variation of only a single 
scalar order parameter - the magnetization. Real interfaces in condensed matter systems are 
much more complex. In a liquid - solid interface for example, in principle, one has to consider the 
variation of an infinity of order parameters corresponding to all the Fourier components of the 
non-uniform solid density p(r). We shall study liquid - solid interface in Chapters 5 and 6. This 
will involve simulations of systems composed of particles, or atoms, with interatomic potentials. 
Before we go on to describe our work on atomic systems, we use this chapter to introduce the 
model interatomic potentials we consider. We then briefly describe the computational techniques 
employed to compute a host of useful thermodynamic, structural and dynamical quantities after 
defining each of them appropriately. Wherever possible, we have compared our results with those 
available in the literature. This constitutes a check on our methods and a verification of the 
codes that we use later to perform simulations of the liquid - solid interface. 



4.2 Models for atomistic systems 

Consider a system containing jV^ atoms. We may divide the potential energy into terms depending 
on the coordinates of individual atoms, pairs, triplets, etc. as : 

^ = Zl^i(''^) + SZ]^2(ri,rj) + ^^ M2(ri,rj-,rfe) + (4.1) 

i i j>i i j>i k>j>i 

J2i^j>i notation indicates a summation over all distinct pairs i and j without counting any 
pair twice; similarly for triplets, etc. The first term X^jMi(ri) represents the effect of an external 
field on the system. The remaining terms represent particle interactions. The second term, U2, 
the pair potential, depends on the magnitude of the pair separation rij — |rj — rj\, so it may be 
written U2{rij). The triplet U3 term is ignored in our thesis The pairwise approximation gives a 
remarkably good description of most atomic systems because the average three-body effects can 



45 



46 



Atomic Systems 



be usually included by defining an "effective" pair potential. To do this, we rewrite Eqn. 
in the form 

V^5^«i(r,)+5^5^t.^^^(r„r,) (4.2) 

i i j>i 

For simplicity, we use the notation u(rij) or u(r) for the pair potential. 



Effective pair potentials 

In our simulations to model atomistic systems we have chosen different pair potentials. We 
discuss them briefly with special emphasis on two-dimensional results. 

• The hard core potential is defined by a pair interaction between two classical particles with 
a non-overlap condition. It is given as 

u"^{r) = oo (r < cr) 

= (a<r) (4.3) 

where a is the diameter of the spheres and r is the distance between the two centers of 
the spheres. A peculiarity of the hard core potential is that it sets a length scale, a, but it 
does not set any energy scale. A configuration of two overlapping spheres cost an infinite 
energy. Since u^^ can only take two values and oo, therefore, for the Boltzmann factor, 
we get, e"'^""^ = e~""*, P = l/ksT. Thus, the temperature scales out trivially. Or, in 
other words, hard objects are athermal; all their structural and thermodynamical properties 
do not depend on temperature so that the density p = N/V (or the packing fraction r], 
which is the ratio of the volume of the spheres and the total accessible volume, V) is the 
only relevant thermodynamical variable. Due to this temperature independence, the hard 
sphere model is the simplest non-trivial model for an interaction. A lot of results exist for 
the hard sphere model (see js^ and references therein) which makes it useful as a reference 
system for systems with more complicated interactions and particle shapes. Moreover, the 
equilibrium thermodynamic properties of the hard sphere model can actually be probed in 
nature by examining suspensions of spherical sterically-stabilized colloidal particles jsof. 

In the two-dimensional case of hard disks the close-packed area fraction is r^c = Txpca"^ /A = 
7i/2^/3 = 0.907... corresponding to a perfect triangular lattice with long-ranged transla- 
tional order. The existence of only two phases, fluid and solid, and a freezing transition from 
one to the other as the packing fraction is changed has been confirmed by extensive simu- 
lations of this model system 0,1611. Hard disks form a key system in our understanding 
of both equilibrium and dynamical aspects of a model liquid-solid interface. 
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• The soft sphere potential is defined as 

M^^(r) = e = ar-" (4.4) 

where z/ is a parameter, often chosen to be an integer, a and e have the dimensions of length 
and energy. The soft-sphere potential becomes progressively "harder" as u is increased. 
Soft sphere potentials contain no attractive part. This retains some of the simplicity of hard 
spheres in the sense that there is again one thermodynamic variable. To understand this, 
note that we can eliminate the explicit temperature dependence by defining a* = (PeY^'^a 
and rewriting Eqn. 4.4 as Pu^^(r) = (a*/ry. The dimensionless excess (over the ideal 
gas) thermodynamic properties then depend only on the single dimensionless state variable 

p* = pa*^ = pa^pe)^/" (4.5) 



Here again we have only two phases but unlike hard systems, soft sphere potential is 
analytic everywhere. In our simulation we use u = 12. For the two dimensional system of 
soft disks (with u = 12), a lot of simulation results exist, including the phase diagram and 
an empirical equation of state [62J. 



The Lennard-Jones 12-6 potential : 

u^J(r) = 4e 



r J \r J 



(4.6) 



This potential has a long-range attractive tail of the form — 1/r^, a negative well of depth e, 
and a steeply rising repulsive wall at distances less than r ~ a. The Lennard-Jones potential 
is extensively used to model atomistic systems and indeed provides a fair description of the 
interaction between pairs of rare-gas atoms 0|. 



4.3 Methods for simulations 

We shall now introduce the various simulation methods that are generally used in simulating 
atomic systems and in particular those that we have used in this thesis to understand the behavior 
of liquid - solid interfaces. We begin by briefly discussing the generic Monte Carlo and molecular 
dynamics methods. Special emphasis has been given to explain the simulation of hard systems 
which require some care during simulations to handle the associated non analyticities. We include 
this section in this thesis mainly for completeness. There exist, of course, excellent textbooks 
0,0,0,161 

wh ch explains these techniques in great detail. 
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4.3.1 The Monte Carlo method 

Consider the N particle system. The classical expression for the thermal average of an observable 
A is given by 

_ /cipV^(r^p-^)exp[-/37i(r^,p^)] 
^ ' /dp^dr^exp[-/3?i(r^,p^)] ^ ' 

where stands for the coordinates of all N particles, and p^ for the corresponding momenta 
and (3 = l/ksT. ?i(r^,p^) is the Hamiltonian of the system. As 7i is a quadratic function of 
the momenta, the integration over momenta can be carried out analytically. The integral over 
particle coordinates is carried out numerically by the Monte Carlo method or, more precisely, the 
Monte Carlo importance-sampling algorithm introduced in 1953 by Metropolis et. al.. 

A basic Monte Carlo algorithm 

In the approach introduced by Metropolis et. al. the following scheme is proposed : 

• Select a particle at random, and calculate its energy U{r^). 

• Give the particle a random displacement; r' = r + A, and calculate its new energy U{v'^) 

• Accept the move from to r'^ with probability 

acc{o — >■ n) = mm(l, exp{-/3[Z^(n) - W(o)]}) (4.8) 

In order to decide whether to accept or reject the trial move, we generate a random number, 
denoted by Ranf, from a uniform distribution in the interval [0,1]. The probability that Ranf 
is less than acc{o n) is equal to acc{o — > n). Therefore, we accept the trial move if 
Ranf < acc{o — > n) and reject it otherwise. This rule guarantees that the probability to 
accept a trial move from o to n is indeed equal to acc{o — > n). 

Using the Monte Carlo method it is particularly easy to to simulate a system of particles 
interacting via hard potentials. The same Metropolis procedure is used, except that, in this 
case, the overlap of two spheres results in an infinite positive energy change and exp{— /9[W(n) — 
U{o)]} — 0. Thus we immediately reject all trial moves involving an overlap since exp{—/5[W(n) — 
U{o)]} < Ranf. Equally all moves that do not involve overlap are immediately accepted. As in 
soft potentials, in the case a move is rejected, the old configuration is recounted in the average. 

4.3.2 The Molecular Dynamics method 

As described in the last section, in a Monte Carlo simulation we compute the average behavior of 
a system in a purely static sense : ensemble average. In most experiments, however, we perform 
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a series of measurements during a certain time interval and then determine tlie average of these 
measurements. The basic idea behind molecular dynamics simulations is precisely that we can 
study the average behavior of system simply by computing the natural time evolution of that 
system numerically. After that we average the quantity of interest over a sufficiently long time. 
The time averaged value of some observable A is therefore given by 

1 /■* 

A{r) = lim - / dt'A{r;t') (4.9) 

t-^oo t Jo 

"Ergodic hypothesis" states that, if we wish to compute the average of a function of the coor- 
dinates and momenta of a many-particle system, we can either compute the quantity by time 
averaging (the "MD" approach) or by ensemble averaging (the "MC" approach). The molecular 
dynamics approach is also useful if we want to look at systems which are driven far from equi- 
librium where it is hard to define an "ensemble" average. We shall use the MD approach for our 
studies of the dynamics of liquids and the liquid solid interface in Chapter 6. 

A basic molecular dynamics algorithm 

• We read in the parameters that specify the conditions of the run (e.g., initial temperature, 
number of particles, density, time step). 

• We initialize the system (i.e., we select initial positions and velocities). 

• We compute the forces on all particles. 

• We integrate Newton's equations of motion. The step and the previous one make up the 
core of the simulation. They are repeated until we have computed the time evolution of 
the system for the desired length of time. To integrate the equations of motion we use the 
"velocity Verlet" algorithm which is implemented in two stages. Firstly, the new positions 
at time t + At are calculated using the equation, 

r(t + At)^r(t)+v(t)At+^At^ (4.10) 

and the velocities at mid step are computed using 

v{t+^At)^v{t) + ^f{t)At (4.11) 

The forces and accelerations at time t + At are then computed, and the velocity move 
completed 

1 1 

v(t + At) = vit + -At) + —fit + At) At (4.12) 

2 2m 
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At this point, the kinetic energy at time t + At is available. The potential energy at this 
time will have been evaluated in the force loop. Note that if we use the kinetic energy 
per particle as a measure of the instantaneous temperature, then we would find that, in 
the canonical ensemble, this temperature fluctuates. To avoid this, we use either simple 
temperature rescaling or proper heat baths j^^l- 

• After completion of the central loop, we compute and print the averages of measured 
quantities, and stop. 

If we have a system of particles interacting via hard potentials then we have to take a different 
approach. Whenever two such particles come close enough to reach a point of discontinuity in 
the potential, then there is a "collision". Therefore, in such a system, the primary aim is to 
calculate the time, collision partners, and all impact parameters, for every collision occurring in 
the system. Unlike soft potentials, hard potential programs evolve on a collision-by-collision basis. 
The general scheme is as follows : 

• locate next collision 

• move all particles forward until collision occurs 

• implement collision dynamics for the colliding pair 

• calculate any properties of interest, ready for averaging, before returning to the first step. 

Although the above scheme is the basis for simulating hard core systems, the introduction of an 
external potential renders the simulation technique difficult to implement. This is because in the 
presence of an external potential the particles will experience a force ftj and the collision equation 
that we need to solve is given by 

\rij{t + ti,)\ = \ri, + ^^ijtij + (1/2)%4| = a (4.13) 

where = rj — Vj and Vj^ = Vj — and a is the hard core diameter. Therefore we now have 
a quartic equation in tij and it is highly non- trivial and time consuming to find the solution. 
Therefore, to simulate a hard disk system in the presence of an external potential, we use multiple 
time step molecular dynamics that we use for soft potentials. However, in the position and velocity 
update scheme, we introduce the additional restriction that if two particles overlap, then they are 
retained in their old positions but their velocities are interchanged. This of course introduces an 
error in the simulation procedure but we choose a sufficiently small time scale and benchmark our 
simulation against standard results for the velocity of sound in hard disk liquid. In our simulations, 
the unit of time, t, is r = ^Jma'^ /ksT , where m(= 1) is the mass of the hard disks. We find 
that a time step of At = 10^"^ conserves the total energy to within 1 in 10^ (at worst) - 10^. 
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4.4 Bulk variables 

From computer simulations, we get information at the microscopic level - atomic and molecular 
positions, velocities, etc. With the help of statistical mechanics we convert this very detailed 
information into macroscopic terms such as pressure, internal energy, etc. In this chapter, we 
focus on the determination of bulk quantities which are averaged over the entire system and are 
meaningful only for homogeneous phases. The analogous space dependent local quantities at 
interfaces are described and discussed in the next chapter. 



4.4.1 Thermodynamic variables 

Bulk thermodynamic variables such as density, pressure, stress etc. can be easily calculated. For 
a system of M particles in a volume V, in a homogeneous phase the average density is given 
simply as p = J\f/V. The pressure may be calculated from the equation ((■) denotes the ensemble 
average) 

VV = N'kBT+{W) (4.14) 
where W is the 'internal virial' given as 

1 ^ 

W = -5^r,-f, (4.15) 

i=l 

in two-dimensions, is the force on some particle i. For pairwise interactions, ^^rj ■ fj = 
X]iX]j>i^«i ■ f« where = ri — rj. We can thus calculate the pressure for various densities. In 
Figs. 14. H and 14. 21 we plot the equation of states for the soft core and Lennard-Jones liquid in two- 
dimensions. In both cases, we perform Monte-Carlo simulations in the constant MAT ensemble 
with J\f = 1200 and T = 1. The potentials are truncated at = 2.5. In our simulations we 
have used reduced units (cr = e = m = 1). The time step of St = 10^^ conserves the energy. 
To compare our simulation results for the soft disk liquid we use the empirical equation of state 
01 given as 

3AV 

= 1.77306p * +2.36241p*2 + 1.798198p *^ -5.648177p *^ 

P 

+78.65712p *^ -197.57241p *^ +212.37417p *^ -79.57456p*^ (4.16) 

For the hard disk system the forces are impulsive and one cannot directly use Eqn. 4.15. We 
however make use of the pair distribution function to calculate the pressure. The pair distribution 
unction is defined in the next section. Thus, the equation of state for the hard disk system is 

Z^(^ = l + 2vg{a) (4.17) 
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Figure 4.1: The equation of state of the soft disk liquid. Symbols are our results, line from |& 
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Figure 4.2: The equation of state of the Lennard-Jones disk liquid. Symbols are our results, line 
from [6£ 



where g{a) is the radial distribution function at r = a, the hard disk diameter. We again perform 
MC simulation in the canonical ensemble with the same number of particles and at the same 
temperature. We average over 10^ configurations to obtain the radial distribution functions for 
various rj's and from them derive the corresponding g{a)'s. We compare our results with the 
semi-empirical equation of state given as j^, 

-1 



1 - 2ri 



2Vc 



Vc 



-7] 



(4.18) 
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Figure 4.3: The equation of state of the hard disk liquid. Symbols are our results, line from 

4.4.2 Structural quantities 

One important function that characterizes the structure of a homogeneous phase is the pair 
distribution function g2{ri, rj), or g2{rij) or simply g{r). It is defined as the probability of finding 
a pair of atoms a distance r apart, relative to the probability expected for a completely random 
distribution at the same density. A form for g{r) useful for computer simulation purpose is given 
as 

9{r) = p-'{J2 E '^(r^)'^(r^- - = E '^(r^- - (4.19) 

In simulations, the delta function is replaced by a function which is non-zero in a small range of 
separations, and a histogram is compiled of all pair separations falling within each such range. 
The pair distribution function is extremely important because we may also express thermodynamic 
quantities like the energy and the pressure in terms of this function. 

One other quantity that is very useful in characterizing the structure of a phase is the structure 
factor S{k) where k is the wavevector. In a simulation with periodic boundaries, k is restricted 
by the periodicity of the system, i.e. with the simulation box. The structure factor describes 
Fourier components of density fluctuations as 

S{k) =Af-\p{k)p{-k)) (4.20) 

where p{k) = Xli^i^^Pl^^ ' ^i)- The structure factor is essentially the Fourier transform of the 
pair distribution function and is directly measurable in X-ray or neutron diffraction experiments 
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Apart from these two quantities we may also define the order parameters that distinguish two 
phases. Here we discuss two such quantities [61] that would be useful to characterize the fluid and 
solid phases in our analysis of the liquid-solid interface of a two-dimensional system in Chapter 
5. 

(i) Bond orientational order parameter : In a periodic crystal, there is only a discrete 
set of directions defined by vectors between nearest neighbour particles, which occupy sites on a 
lattice. These directions are the same throughout the lattice and define a long-range orientational 
order often called the bond orientational order or the bond-angle order. The orientational order 
of a two-dimensional hard disk system can be described by the (global) bond orientational order 
parameter (^g) where as usual (•) denotes ensemble average. The local value of i/jq for a particle 
i located at rj = (x, y) is given by 



where the sum on j is over the Ni neighbours of this particle, and % is the angle between the 
particles i and j and an arbitrary but fixed reference axis. The (global) bond orientational order 
parameter is then defined as 



where is the particle number of the system. For a perfect triangular structure ipe = 1. Thus a 
solid or a hexatic phase will have 'ipQ ^ 0. However for the disordered liquid phase ^Jjq will average 
to zero for A ^ oo. 

(ii) Solid order parameter : In an external field however the bond orientational order may 
be non-zero even in the fluid phase. This is because, we can now have a "modulated liquid" in 
which the local hexagons consisting of the six nearest neighbours of a particle are automatically 
oriented by the external field for example, near the liquid solid interface. Thus ipe is non-zero 
both in the (modulated) liquid and the crystalline phase and cannot be used to study the freezing 
transition or liquid solid interface. The order parameters corresponding to the solid phase are 
the Fourier components of the (nonuniform) density-density correlation (p(rj)p(rj)) calculated 
at the reciprocal lattice points {G}. This (infinite) set of numbers are all zero (for G 7^ 0) in a 
uniform liquid phase and nonzero in a solid. We restrict ourselves to the star consisting of the six 
smallest reciprocal lattice vectors of the two dimensional triangular lattice. In modulated liquid 
phase, the Fourier components corresponding to two out of these six vectors, eg., those in the 
direction perpendicular to the interface, Gi, are nonzero. The other four components of this set 
which are equivalent by symmetry (G2) are zero in the (modulated) liquid and nonzero in the 




(4.21) 
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solid (if there is true long ranged order). Thus, we use the following order parameter : 

N 

(^G,) = (l5^expHGfc-r,,)|) (4.23) 

where = rj — r^. Note that though the order parameter {tpCi) decays to zero with increasing 
system size in the two-dimensional solid - quasi long ranged order - this decay, being weak, does 
not hinder us from distinguishing, in a finite system, a modulated liquid from the solid phase with 
positional order. 

4.4.3 Dynamic quantities 

The important dynamical quantities for our purposes are those which characterize the transfer of 
mass, momentum and energy across a liquid solid interface. In Chapter 6, we shall study the speed 
of sound, the sound absorption coefficient and the thermal conductivity of an inhomogeneous 
system consisting of solid and liquid regions. In this section, we shall show how non-equilibrium 
molecular dynamics simulations may be used to obtain the sound velocity and the absorption 
coefficient for a homogeneous hard disk liquid. Consider a rectangular simulation box of length 
Lx and width Ly. Sound propagation is studied by producing a momentum {vy = Vq) impulse 
over a thin rectangular region in y spanning the width of the cell. This generates a weak, 
acoustic, shock jzd]. The shock rapidly evolves into a Gaussian pulse which propagates non- 
linearly with a speed Cpuise > Co(p), the speed of sound in the liquid which is a function of the 
density p of the hard disk liquid. The deceleration of the pulse follows a pattern which is strongly 




Figure 4.4: Speed of propagation of a momentum pulse in the hard disk liquid as a function of 
time 

reminiscent of the propagation of acoustic shock in water after the collapse of a sonoluminiscent 
bubble. If we fit the time dependent pulse velocity(Fig. 16. 2p to a form Cpuise{t) = co(p) +01/^, 
we obtain the limiting velocity Co(p) with u close to 0.5. We now compare the speed of sound 
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Co(p) with the equation for the speed of sound in a hard disk gas given by j68| 

where for ^(77) we use the familiar equation of state Eqn. 4.18. We find that the speed of sound 
measured in this way compares remarkably well with the speed of sound in the hard disk liquid 
(Fig. 14.511 . The calculated co(p) is of course independent of the initial strength Vq of the shock. 
Viscous dissipation spreads the pulse into a smooth Gaussian with a width which, at late 
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Figure 4.5: Extrapolated limiting pulse speed cq as a function of the packing fraction, 77. The line 
is a plot of the speed of sound obtained from the equation of state Eqn. 4.24 

times, increases linearly with distance s, and the sound absorption [70] coefficient a = A^/4coS 
is also independent of Vq (Fig. 14. 6|) . The g*'^ Fourier component of the velocity therefore decays 
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Figure 4.6: The absorption coefficient showing the expected linear dependence with distance 
of propagation s. The different coloured symbols represent different values of the initial momentum 
shock Vo red : green : etc. 

with s as Vg{s) = t>g(0) exp(— acgg^s). In order to obtain a measurable Gaussian we average all 
quantities over ~ 200 initial velocity configurations chosen from a Maxwellian distribution with 
temperature ksT = 1. 



5 Liquid Solid interfaces : Equilibrium 
properties 



5.1 Introduction and background 

In this chapter we shall use the ideas and methods discussed in the last chapter to study interfaces 
between a liquid and a solid in two dimensions. We shall show that such interfaces may be 
constructed using external chemical potential fields similar to the simple Ising interface studied 
in Chapters 2 and 3. Small wavenumber fluctuations, known as crystallization waves, which are 
dominant in field free liquid solid interfaces are, of course suppressed for reasons similar to the 
Ising case viz. the explicit breaking of translational invariance of the system perpendicular to the 
interface. Our interfaces therefore remain planar for the parameters used in our study. However, 
as we shall see shortly, the planar liquid solid interface is not inert. Particles are transferred 
across the interface in new and interesting ways. We believe that some of our predictions may 
be directly checked for liquid solid interfaces in atomic, as well as, colloidal systems where the 
chemical potential field may be provided either by a laser trap or by adsorbing colloidal particles 
on a patterned substrate. In the next chapter we shall discover how these interfacial fluctuations 
influence transport properties of the interface contributing to the resistance offered to the transfer 
of momentum impulses and of heat. 

A liquid is homogeneous and isotropic. The local density p(r) is independent of the spatial 
coordinate r and has Fourier components which are nonzero only for the q = mode. This 
translational invariance is partially broken in the solid which has Fourier components pc which 
are nonzero for all wavevectors G which are members of the reciprocal lattice vector set of the 
solid. The solid can therefore be described completely using the countably infinite number of 
Fourier components pc- Across a liquid solid interface all the order parameters are expected to 
go from zero to their value in a solid (close to unity) within a region of width ^, the interfacial 
thickness. While, the structure of real liquid-solid interfaces is difficult to determine experimentally 
in molecular detail, the interface is expected to be sharp with the order parameters building up 
within a few atomic layers. This is simply a reflection of the fact that the freezing transition 
is sharply first order in three dimensions j71| with a finite correlation length. Also most liquid 
solid interfaces are not planar. The finiteness of heat and mass transport coefficients at the 
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interface generates instabilities, first analyzed by Mullins and Sekerka p, |72|, which grow into 
non-compact dendritic structures typically observed in most real situations. There has, however, 
been a number of theoretical and computer simulation studies starting from the pioneering work 
of Haymet and Oxtoby iZ4i]. Most of these studies refer to a planar interface. 

In the next section we shall introduce the local analogs of the thermodynamic and structural 
quantities introduced in the earlier chapter. All of these quantities will have different values on 
either side of a liquid solid interface and will be used to define the interface. Next we shall show 
how non uniform chemical potential fields may be used to generate liquid solid interfaces and 
shall try to make contact with experiments by calculating approximately the parameters required 
for constructing a suitable laser trap. We shall plot local thermodynamic and orderparameters 
to characterize the interface. Our main result of this chapter, namely the existence of complete 
layer transfers will then be presented for three different systems viz. the hard disk, soft disk 
and Lennard-Jones system. We shall attempt to understand our result from thermodynamic 
arguments, for the simplest of the three systems, namely hard disks. 



5.2 The liquid solid interface 

in this section, we explore the possibility of creating a patterned sequence of confined solid and 
liquid regions using an external, space-dependent, chemical potential field </)(r). Consider a two 
dimensional system (see Fig 15. Ij) of N atoms of average density (packing fraction) i] = -kN/AA 
within a rectangular cell of size A = ^ Ly where the central region S of area As = x Lg 
is occupied by A^^ atoms arranged as a crystalline solid of density rfs > rj, while the rest of the 
cell is filled with liquid of density rji < r]. The difference in density is produced by an external 




Figure 5.1: A cartoon diagram of the system 



field 0(r) = — /i for r G 5; increasing sharply but smoothly to zero elsewhere with a hyperbolic 
tangent profile of width S^f^. How may 0(r) be realized in practice? In model solids like colloids 
|75|, one may use a surface template to create a static pattern In real systems, as well as 
colloids jzyi, one may be able to use laser traps j3| or non -uniform electric or magnetic fields. 
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Usual laser traps for alkali metals or rare gas atoms are in the range of lOmK for which a power 
of about lOOmW is required. 

For specificity and simplicity we first choose the atoms to interact with a hard disk potential 
lygj l although we show towards the end of this chapter that qualitative results for more realistic 
potentials, e.g., soft sphere or Lennard-Jones are the same. We have chosen = a/ A, where 
a is the hard disk diameter and sets the scale of length. The energy scale for this system is set 
by ksT where ks is the Boltzmann constant and T the temperature. In our simulation we set 
a = ksT = 1. Most of the qualitative phenomena discussed subsequently are independent of 
the exact nature of the interactions. 

The full configuration dependent Hamiltonian is 7i = y^, . + J^i'Pi^i)- \^3ve carried 
out extensive MC simulations with usual Metropolis moves (64l, periodic boundary conditions 
in both directions and in the constant number, area and temperature ensemble to obtain the 
equilibrium behaviour of this system for different /i at fixed r]. N = 1200 particles occupy an 
area A = 22.78 x 59.18 with the solid occupying the central third of the cell of size Ls = 19.73 
The initial configuration is chosen to be a liquid with i] = .699; close to but slightly lower than 
the freezing jz^ density rjj = .706. On equilibration, S contains a solid with the closest -packed 
planes parallel to the solid -liquid interfaces which lie, at all times, along the lines where </)(?/) 0. 
The equilibration time is large and many (~ 10'') Monte Carlo steps (MCS) are discarded before 
results shown in Figs. I5.2l and 15.31 are obtained. 



5.3 Local averages 

Since in this chapter we study a liquid solid interface, thermodynamic and structural quantities 
which are averaged over the entire system makes no sense. We have to, therefore, construct local 
analogs of these quantities as discussed in Chapter 4. To this purpose, we divide the simulation 
box into thin rectangular strips with the longer dimension parallel to the interface. Each of these 
strips are indexed by an integer i which runs from 1 to M, the total number of strips. We shall 
see that the external potential induces a solid which has a layered structure. Therefore, it is 
important that while breaking the system up into strips of equal width, the edge of the strips in 
the solid region always falls between two layers. The index i multiplied by the width of the strip 
gives the distance in the direction perpendicular to the interface. 

The thermodynamic quantities may always be defined as intensive variables containing sums 
over all the particles of the system divided by the total number of particles. In this case the sums 
run only over the particles in a strip. Thus for example. 

Number of particles in strip i 
Strip width 

where pi is the local density in strip i, and similarly for other quantities. 
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0.76 




Figure 5.2: The density profile ri{y) coarse grained over strips of width a (averages taken over 10^ 
MC configurations each separated by 10^ MCS), varies from r]£ to rjg as we move into S. The lines 
show the predictions from a simple thermodynamic theory presented in Section 5.5.1 

Structural quantities are typically sums over pairs of particles. Here, we take one of the 
particles within the strip while its companion lies either in the same strip or in neighbouring strips. 
The orientational order parameter is calculated in each strip and is averaged over equilibrium 
configurations. To calculate the local solid order however, we choose area within the solid region 
and count the number of particles in this area. We then find out the solid order parameter for 
these particles as a function of the length perpendicular to the interface. We use the same trick 
for the liquid region and finally merge the two at the interface to get the variation of the solid 
order parameter for the entire system in the direction perpendicular to the interface. Just as for 
the orientational order parameter, averaging is done over equilibrium configurations. 

In all these definitions, we have of course assumed that there is no systematic variation of any 
of the quantities in the direction parallel to the interface, in accordance with the symmetry of 
the problem. 

5.4 Order parameter profiles and structure factor 

In this section we present all results for the local thermodynamic and structural quantities calcu- 
lated as above. The density r]{y) coarse grained over strips of width ~ a varies from its value r]i 
in the liquid to r]s as we move into (and away from) S (Fig. I5.2jl . Averages taken over 10^ MC 
configurations each separated by 10'^ MC steps. The trap depth /i = 6, supports an equilibrium 
solid of density rjs = .753 in contact with a fluid of density r]i = .672. The horizontal lines 
are predictions of a simple free-volume based theory (discussed in later sections) for 77^ and r]i. 
A superposition of atomic positions shows a static, flat, liquid solid interface with the solid like 
order gradually vanishing into the liquid (Fig. 15. 3p . We have thus created a thin nano-sized 
crystal which is about 21 atomic layers wide (for a trap depth, /i = 6) and is flanked on either 
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side by liquid separated by two liquid solid interfaces. 




Figure 5.3: Superposition of particle configurations from the MC run showing a solid like order (red 
: high rj) gradually vanishing into the fluid (yellow : low rj) across a well defined solid-fluid interface. 

The bond orientational order parameter {ipe{y)) coarse grained over strips of width a (averages 
taken over lO'^ MC configurations each separated by 10^ MCS) shows a sharp rise from zero to a 
value close to one, as we move into the region S. This indicates that the particles in this region 
are arranged approximately hexagonally. However, this does not necessarily justify the phase to 
be a solid. Therefore, the need to calculate the solid order parameter. The solid order parameter 
(i/jg) in the direction G2 and the other equivalent to it by symmetry is nonzero in the region S, 
(Fig. I5.5jl indicating without doubt the nucleation of a solid phase. Note that (^g) shows a larger 
interfacial region than that obtained from (^^02)- This is because the liquid near the liquid solid 
interface is highly orientationally ordered. 

The two-dimensional structure factor shows sharp Bragg peaks for the solid region and isotropic 
pattern for the liquid (Fig. 15.611 . Although the structure factor for the solid shows well defined 




Figure 5.4: Bond orientational order parameter across the liquid solid interface for a 21 layered 
solid surrounded by liquid on both sides. 
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Figure 5.5: (a)The reciprocal lattice vectors Gi and G2, and the rectangular unit cell, (b) Solid 
order parameter corresponding to G2, is nonzero in the region S, indicating a solid 



peaks with six-fold symmetry of the triangular lattice, the quasi-onedimensional nature of our 
system implies that true solid -like order may not be present. However, we do not pursue this 
question in detail in this thesis. In the interfacial region the peaks are seen to be diffuse indicating 
once more that the interfacial region is highly orientationally ordered. 
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Figure 5.6: Two-dimensional structure factor for the liquid, solid and interfacial regions respectively. 



5.5 The layering transition 

We now calculate the difference in densities between the solid and liquid regions Ar] = {rjs — rfi) 
as a function of the strength of the external field /i. While Arf/rj increases with increasing fi 
as expected, the smooth increase is punctuated by a sharp jump (Fig. 15. 7j) . An examination of 
the particle configuration shows that the jump occurs when an extra close -packed layer enters 
S increasing the number of solid layers by one. For the parameters in our simulation, the jump 



5.5 The layering transition 



63 







C 


^^^^^ 


- 






B 











2 4 ,, 6 8 10 



Figure 5.7: Plot of the equilibrium fractional density change Ary/ry as a function of ii (points (MC 
data), thick solid line (approximate theory)), showing discontinuous jump at /i « 8. 



occurs at /i ^ 8 with the number of layers increasing from 21 to 22. The value of at the 
jump is a strong function of Lg. The solid structure is seen to be a defect free triangular lattice 
with a small rectangular distortion SdirisjLs) [8q|. We have examined the variation of Ari{jj,)/ri 
by cycling fi adiabatically around the region of the jump. This yields a prominent hysteresis 
loop as shown in Fig. I5.8l which indicates that 'surface' steps (dislocation pairs) nucleated in the 
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.11 
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10 



Figure 5.8: Cycle-averaged hysteresis loop as fi is cycled at the rate of .2 per 10^ MCS. The central 
jagged line is the result of the initial cycle when a single dislocation pair was present in the central solid 
region. 

course of adding (or subtracting) a solid layer, have a vanishingly short lifetime. Consistent with 
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this we find that the jump in Ari{fi)/ri vanishes when the system is minimised at each fi with a 
constraint that the solid contains a single dislocation pair (Fig. 15. 8p . Interestingly, a dislocation 
pair forced initially into the bulk, rises to the solid-fluid interface due to a gain in strain energy 



j8ll |. where they form surface indentations flanked by kink-antikink pairs. This costs energy due to 
the confining potential, as a result, the kink-antikink pair gets quickly annealed by incorporating 
particles from the adjacent fluid. The jump in Ar]{fi)/r] is also seen to decrease with increasing 



5.5.1 The thermodynamic theory 

The qualitative features of these results may be obtained by a simple thermodynamic theory 
(Fig. 15. 7 j) with harmonic distortions of the solid, ignoring contributions from spatial variations 
of the density, thermodynamic theory. We first write down the free energy densities of the bulk 
fluid and solid phases. 

Free energy of the solid 

For a solid which is only a few atomic layers wide, it is impossible to separate surface and bulk 
contributions to the free energy. We therefore assume that the dominant effect of confinement of 
the solid is the introduction of an uniform strain to be calculated from a reference triangular 
lattice with the same number of atomic layers [80.]. Therefore, we have two terms in the free 
energy of the bulk solid. 

(i)Pree volume part : In a given solid lattice, the particles have a bulk mean particle distance 
tto which is the distance between nearest neighbours of the lattice. We consider the particles 
to be confined in independent Wigner-Seitz (or Voronoi) cells of the solid. The cell has a 
volume Vb = gbO-o where gb is a geometrical prefactor that depends on the lattice type and 
on dimensionality D. Each center-of-mass coordinate of the hard disks can move within a free 
volume of 



without touching the neighbouring disks. Hence, one obtains a lower bound for the bulk partition 
function, Q > (V;,//A^)-^ which provides an upper bound to the bulk free energy density (free 
volume contribution) 



The upper bound becomes asymptotically exact for close packing. In our calculation for the 





Vbf = Qbia-Q - 



(5.2) 




(5.3) 



two-dimensional hard disk solid we calculate gb = 1/(0.5231(t)^ and aQ = a 
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(ii)Elastic contribution : In order that the solid channel accommodates rii layers of a homo- 
geneous triangular lattice with lattice parameter of hard disks of diameter cr, we need 



V3 

— {ni-l)ao + a (5.4) 



Defining 



X(r/s, ^.) = 1 + (5.5) 

Eqn. 15.41 then implies x = integer = rii and violation of Eqn. 15.41 implies a rectangular strain 
away from the reference triangular lattice of rii layers. The lattice parameters of a centered 
rectangular lattice (CR) unit cell are and ay (Fig. 15.511 . In general, for a CR lattice with given 
Ls we have ay = 2(Ls — cr)/ [rii — 1) and, = 2/ pay. The normal strain Ed = Sxx — £yy is then, 

111 — 1 X ~ ^ 
Sd = -, (5.6) 

X - 1 rii-l 

where the number of layers rii is the nearest integer to x so that Ed has a discontinuity at half- 
integral values of x- For large Lg, this discontinuity and Ed itself vanishes as l/Lg for all rjs. We 
therefore have an elastic contribution of {l/2)KA£d in the free energy of the solid. 

Incorporating these two factors we may write down the free energy density of the solid phase 

as 

fs = fA + \k^eI (5.7) 
where /Ca is the Young's modulus of an undistorted triangular lattice. 

Free energy of the liquid 

The free energy density of the liquid bulk phase may be simply written as 

fi = r dv[^^^^ + hd (5.8) 

where the ideal gas Helmholtz free energy per particle fid = ln(p) — 1 and we use the semiemperical 
equation of state for the hard disk liquid [gQ^ given by Eqn. 14.181 

Free energy of the system 

We now write down the total free energy density of the system (fluid + solid regions) using the 
free energy density expressions for solid and liquid bulk phases as 



f = x [fsiVs, Ls) - %/i/vr] + (1 - x)fi{r]t,) 



(5.9) 
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We then minimise this free energy density with the constraint that the average density is fixed, 
T] = xr]s + (1 — x)r]£, where x is the area fraction occupied by S. The result of this calculation 
is shown in Fig. I5.7l where it is seen to reproduce the jump in A?7(yu)/?7. 

Why does the solid incorporate layers of atoms from the liquid ? This question may be answered 
elegantly if one calculates the uniaxial stress in the solid region as a function of the depth of the 
strain, e^. The stress may in fact be obtained in a straight forward fashion from the expression 
of the free energy. Differentiating the free energy of the solid w.r.t. Ed we obtain 

^-t (5.10) 

When 7d is plotted versus the uniaxial strain Ed, we observe that the solid is not stress free for any 
arbitrary value of a combination of /i and Lg. In fact, for our parameters, initially the 21 layered 
solid is under tension. We follow the variation of the tensile stress with the strain as fx is increased 
from the points A — D \f\ the Fig. 15.91 The state of stress in the solid jumps discontinuously 




Figure 5.9: A plot of the tensile stress 7^ against strain Ed- The arrows show the behaviour of 
these quantities as fi is increased from the points marked A - D. 

from tensile to compressive from B^C due to an increase in the number of solid layers by 
one accomplished by incorporating particles from the fluid. This transition is reversible and the 
system relaxes from a state of compression to tension by ejecting this layer as n is decreased. As 
fi is increased, the tension increases till it reaches about —2.45 when the corresponding strain 
is about —0.038. At this point a layer enters the solid region and the stress and strain switches 
from tensile to compressive. Further increase in yU will now decrease the stress and drive the solid 
to a state of zero stress. Thus, the layering transition from 21 to 22 layers as observed by us is 
a mechanism for relieving stress. 

Solids subject to large uniaxial deformations, relieve stress either b y th e ^generation and mobility 
of dislocations and/or by the nucleation and growth of cracks (la, |S2|- What is the nature 
of stress relaxation when conditions are arranged such that these conventional mechanisms are 
suppressed ? Nano-indentation experiments jssi, Is^ show that if small system size prevents the 
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Figure 5.10: Plot of the equilibrium fractional density change Ap/p as a function of fi (points 
(MC data) for soft core; line denotes approximate fit 



generation of dislocations [S(^, solids respond to tensile forces by shedding atoms from the surface 
layer. What we have shown here using extensive computer simulations and theory, is that a small 
crystal trapped within a potential well ||85| and in contact with its own fluid, a situation easily 
realised using optical traps, responds to large compressive stresses by a novel mechanism — the 
transfer of complete lattice layers across the fluid solid interface. In Chapter 6, we study the 
importance of this "layering" transition by looking at momentum and heat transport across the 
liquid solid interface and the role this transition plays on transport properties. 

5.5.2 Layering in other potentials 

If the layering transition observed in the hard disk system is actually a new nechanism for relieving 
stress in a thin crystal, it should be independent of the details of the potential. Our main results 
trivially extend to particles interacting with any form of repulsive potential, or even when the 
interactions are augmented by a short range attraction, provided we choose fi deeper than the 
depth of the attractive potential. In this section we show explicitly that the layering transition 
is present in the soft core and Lennard-Jones system. Once again we perform MC simulations in 
the constant MAT ensemble with periodic boundary conditions and with the external chemical 
potential /i. The relevant parameters corresponding to these potentials, namely e and a in Eqn. 
4.4 and Eqn. 4.5, set the energy scale and the length scale respectively. In our simulation 
e = cr = 1 and J\f = 1200 particles occupy an area ^ = 24 x 60 with the solid occupying 
the central third of the cell of size Lg = 20. The average density of the system is therefore 
p = 0.833 to be compared with the freezing density of p ~ 1.0 and p ~ 1.0 for the soft core and 
Lennard-Jones systems at this temperature. For the soft core potential, the fi value at the jump 
in density (Fig. IS.lOp is even quantitatively comparable to the corresponding hard core system. 
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Figure 5.11: Plot of the equilibrium fractional density change Ap/p as a function of fi (points 
(MC data) for Lennard-Jones system; line denotes approximate fit 

In the LJ system the jump happens earlier as expected due to the cohesive energy, Ccoh of the LJ 
atoms (Fig. 15.111) . We estimate, /i' = ^hd — Gcoh, where ^hd is the jump potential for the hard 
disk. Thus Ccoh ~ —3.5. The shift of the layering transition may be used to obtain the cohesive 
energy of the trapped solid. 




6 Liquid solid interfaces : Dynamical 
properties 



6.1 Introduction 

In the last chapter, we showed how one might produce a liquid solid interface using a non- 
uniform external potential. We characterized this interface using a variety of thermodynamic 
and structural quantities, which were measured as a function of the perpendicular distance from 
the interface. Finally, we showed that as a function of the depth of the potential well, the 
trapped solid undergoes, what we called, "layering" transitions which involved the addition (or 
removal) of an entire layer of solid from or into the surrounding liquid through the interface. The 
layering transition was accompanied by a sharp jump in the density of the solid. We showed that 
this layering transition is a novel mechanism by which a stressed nano-solid constrained by an 
external potential can respond plastically to large stresses without nucleating dislocations. Finally, 
we established that this phenomenon is general and is independent of the particular interatomic 
potential used. 

In this final chapter of the thesis we shall explore how mass, momentum and energy is trans- 
ported across the liquid solid interface and especially the role of the layering transition on the 
transport coefficients. We shall show in this chapter that (1) fluctuations associated with this 
transitions are of a special kind always involving the transfer of complete layers of solid (2) these 
fluctuations offer resistance to the transfer of momentum and energy through the interface (3) 
the resistance is maximum when the energy matches that required to raise a complete lattice 
layer from within the potential well into the surrounding liquid. In the next section we shall begin 
by studying the stability of surface kinks at the liquid solid interface in the hard disk system. 
We shall then study the response of the interface to weak acoustic shocks. In the last section 
we shall use non-equilibrium molecular dynamics to study heat transport through the liquid solid 
interface in soft disks and obtain the contact or Kapitza resistance of the interface as a function 
of the depth of the potential well. 

In our investigations, to be reported in this chapter we shall exclusively use molecular dynamics 
simulations since we shall be interested in dynamical quantities. 
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6.2 Kinetics of layering 

The large hysteresis loops associated with the layering transition obtained in Chapter 5 makes it 
clear that the kinetics of this transition is slow. To study the lifetime of the kink-antikink pairs 
(surface step), we resort to a MD simulation. 

Starting with an equilibrium configuration for the hard disk system, taken from our Monte 
Carlo runs as discussed in the last chapter (Section 5.2) at /i = 9.6 corresponding to a 22-layer 
solid, we create a unit surface step of length / by displacing a few interfacial atoms from the 
solid region into the liquid and 'quench' across the transition to fi = 4.8, where a 21-layer solid 
is stable. The rest of the parameters are kept identical to those given in chapter 5. We observe 
that the fluctuation thus created rapidly relaxes back and the surface step vanishes as the atoms 
are pulled back into the solid. We illustrate this by plotting the number of hard disks within the 
solid region as a function of the MD time steps (Fig. 6.1). As soon as a step is created, the 




t 

Figure 6.1: Plot of the number of particles in the solid region {Ms) as a function of time t clearly 
shows that the displaced particles are pushed back into the region. 

line of atoms in the portion of the solid thus exposed bend to fill up the gap created between 
the atomic layers and the edge of the potential trap. This generates considerable local elastic 
stress. Also, the liquid layer lying immediately adjacent to the solid has a lot of orientational 
and solid-like order. For short times it responds elastically to the presence of an increased local 
density of atoms. The combined elastic response therefore pushes the displaced atoms back into 
the solid region thereby annihilating the step. Indeed, a free energy audit involving a bulk free 
energy gain AF ~ l/Ls, going from a 22 to 21 layered solid, and an elastic energy cost ~ log(/) 
for creating a step of size /, reveals that a surface step is stable only if / > /* ~ l/Ls- For small 
Lg, the critical size /* may therefore exceed L^, the total length of the interface. Of course, if 
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the step spans the entire length of the interface, there is no bending of the atomic lines and there 
is no elastic energy cost. This explains the slow kinetics since the system has to wait till a rare 
random fluctuation, which displaces all the atoms in a solid layer across the interface coherently, 
is required for the layer transition to occur. Although we have explicitly demonstrated this for 
the hard disk system, we believe that similar considerations should be appropriate for the soft 
disk and Lennard Jones systems too. 

The slow kinetics of the layering transition may have an impact on the transfer of momentum 
across the liquid solid interface in the form of regular sound waves or acoustic shocks. The large 
effective compressibility of the solid at the layering transition as evidenced by the jump in the 
density as the chemical potential is increased by an infinitesimal amount (Fig. I5.7jl should reduce 
the velocity of sound considerably. The propagation and scatterin g of sound in an inhomogeneous 
region with coexisting phases has been studied extensively in the past. The 

transfer of mass between coexisting phases at inter-phase boundaries is known to slow down 
and dissipate sound waves traveling through the system. Our system has an artificially created 
inhomogeneity, which should have a similar effect on its acoustic properties. 

Further, the mechanism of stress relaxation of a thin {Lg small) solid via the transfer of an 
entire layer of atoms may be exploited for a variety of practical applications. Provided we can 
eject this layer of atoms deep into the adjoining fluid and enhance its lifetime we may be able 
to use the ejected layer of atoms to create monolayer atomic films or coatings. Highly stressed 
mono-atomic layers tend to disintegrate or curl up [89J as they separate off from the parent 
crystal. It may be possible to bypass this eventuality, if the time scale of separation is made much 
smaller than the lifetime of the layer. Can acoustic spallation [90] be used to cleave atomic layers 
from a metastable, stressed nanocrystal? 

In the next section, therefore, we study the response of the liquid solid interface in our system 
of hard disks to acoustic shocks with a view to studying the effect of the layering transition on 
acoustic shock propagation and dissipation. 

6.3 Effect of shock wave 

imagine, therefore, sending in a sharp laser (or ultrasonic) pulse, producing a momentum impulse 
{vy{t = 0) = Vo) over a thin region in y spanning the length of the simulation cell, which 
results in a weak acoustic shock j^B] (corresponding to a laser power ^ 10^ mW and a pulse 
duration Ips for a typical atomic system). The initial momentum pulse travels through the solid 
and emerges at the far end (Fig. I6.2jl as a broadened Gaussian whose width. A, is a measure of 
absorption of the acoustic energy of the pulse due to combined dissipation in the liquid, the solid 
and at the interfaces |j70|, ||86|, i9ij . For large enough pulse strengths Vq, this is accompanied by 
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Figure 6.2: (a)-(c) Plot of the absolute value of the momentum for molecular dynamics 

times t = .0007 (a), .2828 (b) and 2.8284 (c). The green lines show the position of the solid -fluid 
interfaces. The fit to a Gaussian (blue line) is also shown in (c). Curves such as in (a)-(c) are obtained by 
averaging over 100 — 300 separate runs using different realizations of the initial momentum distribution. 
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a coherent ejection of the (single) outer layer of atoms into the fluid. Note that in our MD 
simulations, to reduce interference from the reflected pulse through periodic boundary conditions, 
we increase the fluid regions on either side, so that for the MD calculations we have a cell of size 
22.78 X 186.98 comprising 3600 particles. In Fig. I6.2l the initial momentum pulse with strength 
Vo = 6. is given within a narrow strip of size ~ a, just to the left of the solid region and the 
curves are fitted to a Gaussian (and the width extracted) when the maximum of the pulse 
reaches a fixed distance of 44.1 from the source. A reflected pulse can also be seen. 

When a shock wave, which propagates through a conventional solid emerges from the free 
surface, the compressed material expands to zero pressure jo^]. The unloading (rarefaction) 
wave travels backwards into the material with the speed of sound. The response of the solid 
now depends on the specific nature of the shock front. For a shock wave with an approximately 
Gaussian profile as in our case, significant negative pressures can develop at the interface where 
the shock emerges due to the interaction of the forward and the reflected waves and a portion 
of the solid may split off by a process known as "spallation". Spallation in bulk solids like steel 




Figure 6.3: Plot of the squared width of the momentum pulse after it emerges from the solid 
as a function of Vq for /x = 4.8(0) and 9.6(n). The solid line is a guide to the eye. The peak in 
A^(Vo) so produced is more prominent for the metastable 22 layered solid /i = 4.8 than for the stable 
[fj, = 9.6) system showing a more coherent momentum transfer in the former case. 
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needs acoustic pressures in excess of 10^ N/cm^ ^] usually available only during impulsive 
loading conditions; the ejected layer is a "chunk" of the surface. In contrast, the pressures 
generated by the shock wave in our system causing coherent nano spallation involves surface 
stresses of the order kBT/a"^ ^ 10^^ N/cm^, and are therefore much lower. This difference 
comes about because unlike a bulk system, a strained nanocrystal on the verge of a transition 
from a metastable n + 1 to a n layered state readily absorbs kinetic energy from the pulse. The 
fact that surface indentations are unstable (Fig. 16. 4p unless of a size comparable to the length 
of the crystal, L^, ensures that a full atomic layer is evicted almost always, leading to coherent 
absorption of the pulse energy. The coherence of this absorption mechanism is markedly evident 
in a plot of against Vq which shows a sharp peak (Fig. 16. 3p . Among the two systems studied 
by us, viz., a metastable (/i = 4.8) and a stable (/i = 9.6) 22 layered solid, the former shows a 
sharper resonance. Note that the absorption of momentum is largest when the available kinetic 
energy of the pulse exactly matches the potential energy required to eject a layer. To elucidate 
this fact further, we plot the configurations of the metastable system (/i = 4.8) as the pulse 
travels through the system, for two different pulse strengths, Vq = 2 and Vq = 6 (Fig. 16. 4|) . The 
weaker momentum pulse (Vq = 2) initially ejects a few atoms of the interfacial crystalline layer of 
the metastable 22 layered solid. However, the resulting large non-uniform elastic strain evidenced 
by the bending of lattice layers causes these atoms to be subsequently pulled back into the solid. 
This effect is the same as that seen in Section 6.2. Only a stronger pulse, Vq = 6, capable of 
ejecting a complete lattice layer succeeds in reducing the number of solid layers by one leading 
to overall lower elastic energy. 

The eviction of the atomic layer is therefore assisted by the strain induced interlayer transition 
and metastability of the 22 layered solid discussed above. Spallation is also facilitated if the atomic 
interactions are anisotropic so that attraction within layers is stronger than between layers (eg. 
graphite and layered oxides [89]), for our model, purely repulsive, hard disk solid, an effective, 
intralayer attractive potential of mean force is induced by the external potential [8Q||. 



The spallated solid layer emerges from the solid surface into the fluid, and travels a distance 
close to the mean free path; whereupon it disintegrates due to viscous dissipation (Fig. 16. 5p . 
To obtain the life time of the spallated layer we calculate the time development of the Fourier 
component of the local density correlation pG(y,t) obtained by averaging, at each time slice t, 
the sum J2j=iN^^P(~^^' i^i ~ '^i)) ^" P^^icle positions ri within a strip of width ~ o 
centered about y and spanning the system in x. The wavenumber G = (27r/(i)n where d = .92 
is the distance between crystal lines in the direction fi normal to the fluid-solid interface. The 
solid (central region with pQ,{y,t) ^ 0) ejects a layer (shown by an arrow in Fig. 16.511 which 
subsequently dissolves in the fluid. The lifetime of the layer is around 2-3 time units (r) which 
translates to a few ps for typical atomic systems. The lifetime increases with decreasing viscosity 
of the surrounding fluid. Using the Enskog approximation j92| to the hard disk viscosity, we can 
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Figure 6.4: (a) Configuration snapshot from a portion of our MD cell showing hard disk atoms 
(green circles) at the solid (bottom)- liquid interface (yellow line) as a weak momentum pulse (Vq = 2.) 
emerges into the liquid, at three different times, for jj, = 4.8. The pulse initially ejects a few atoms 
from the solid but are subsequently pulled back due to the large elastic strain cost in bending, of the 
interfacial crystalline layer (red circles) of a metastable (b) The same for a stronger momentum pulse, 
Vb = 6. This time the pulse strength is sufficient to eject the layer. 
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115 



Figure 6.5: A plot of the time development pQ,{y,t) The solid ejects a layer (shown by an arrow) 
which subsequently dissolves in the fluid. The curves from bottom to top correspond to time slices at 
intervals of At = .07 starting from t = 1.06 (bottom). We have shifted each curve upward by .03t/At 
for clarity. 



calculate the bulk viscosity for a hard disk fluid to be 



001 



Qe = —CooV 9{(^) 

TT 



(6.1) 



where Coo is a constant and g{a) is the pair-correlation function at contact which can be calculated 
using Eqn (4.17) and Eqn. (4.18). For a system of hard disks with m = a = /3 = 1, we can 
calculate Coo = l/2v/7r- Thus, (e oc rf and we estimate that by lowering the fluid density 
one may increase the lifetime of the layer considerably. The lifetime enhancement is even greater 
if the fluid in contact is a low density gas (when the interparticle potential has an attractive part 

H). 
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Figure 6.6: A plot of the total number of particles Ng within the solid region [fi = 4.8) as a 
function of time for Vq = 1. (top) and 6.. Note oscillations in Ng] only the stronger pulse changes the 
number of solid layers from 22 to 21. 

6.4 Effective liquid theory 



The absorption line -shape may be understood within a phenomenological "effective liquid" ap- 
proximation. The extra absorption producing the prominent peak in Fig. 16.31 is due to the loss of 
a whole layer from the solid into the liquid (Fig. I6.4jl . For small Vq, the confined solid responds 
by centre of mass fluctuations {q ^ phonons) shown by oscillations of A^s with time (Fig. 16. 6j) . 
Scattering from this and other sources 0, m 0, H 01 constitute a background which we 
ignore, as a first approximation, for simplicity. Within our approximation, the momentum loss 
at the interface is modelled as regular dissipation within a liquid strip of (fictitious) width ^ (see 
Fig. I6.7jl . The expected momentum transfer at the interface Ug = UqX the probability that the 



11 



local 



Figure 6.7: A cartoon diagram showing the momentum transfer as assumed in our phenomenolog- 
ical theory 
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momentum Ho required to eject the layer, exists. If a local "temperature" Tiocai measures the de- 
gree of (de)coherence of the momentum transfer, then lie = (l/2)noerfc[(no — Vq)/ y/2kBTiocai] 
and ^ may be extracted from Vq — He = Vq exp(— au;^^). Substituting for Hg we obtain the extra 
absorption due to the interface, 

Vo 



= A.acl^ = -alog[l - ^erfc{^^= 



}] 



(6.2) 



We use a, Ho and Tiocai as fitting parameters. In Fig. 16.31 we show a fit to Eq. 16.21 of our MD 
data and observe that it reproduces all the features remarkably well. The larger error-bars near 
the peak in A^(Vo) reflects the difficulty of fitting a Gaussian to the transmitted pulse when 
dissipation is large. Indeed, in this region the pulse shape is systematically distorted away from 
Gaussian due to effects beyond the scope of our simple theory. Large fluctuations in A^.. (Fig. l6.6|l 



lead to expected 



rects Devon a tne s 
Zq|,H0,H| an 



d detectable decrease in average pulse speed (Fig. I6.8jl . 




Figure 6.8: Average pulse velocity Cpuise{Vo) for = 4.8; note the dip in Cpuise where absorption 
is strongest. 



6.5 Heat transport across the liquid - solid interface 

in the last section we studied the effect of the layering transition on the transport of momentum 
across the liquid solid interface. In this section we shall focus on the transport of energy. The 
kinetic energy flux jg is related to the temperature gradient VT through the famous Fourier's 
law, viz. 

ie = -AVr (6.3) 

where, the transport coefficient A is the thermal conductivity. The transport of heat through 
small and low dimensional system has enormous significance in the context of designing useful 
na no-structures [91j. Heat transport across a model liquid solid interface has been studied in 
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three dimensions with the interatomic potential being Lennard-Jones [95j. They have shown that 
the Kapitza resistance can reach appreciable values when the liquid does not wet the solid. In 
two dimensions, strictly speaking, Fourier's law needs to be modified by correction terms of the 
order of the logarithm of the system size [96] . For a small system, however, we shall ignore these 
corrections and assume it to be valid. 

The specific context in which we study thermal properties of the liquid solid interface is the 
same as that we have used for our studies of the acoustic properties in the last section. Again, 
a solid region is created within a liquid using an external chemical potential trap. The overall 
dimensions of the system and the size of the solid region are also identical. The inter-atomic 
potential however is chosen to be soft sphere with the potential given by Eqn. 4.4. The smooth 
potential allows us to use standard MD simulations with a velocity Verlet algorithm. The time 
step of St = 10^^ in our MD ensures that the total energy is conserved (in equilibrium) to 
within 10"'^. Periodic boundary conditions are applied in the x-direction. A temperature gradient 
is generated throughout the simulation cell by coupling the two sides of the system in the y- 
direction to two reservoirs at different temperatures, which drives the system out of equilibrium. 
Thus we have a reservoir at y = at a temperature Ti and another at y = Ly at temperature 
T2. The coupling to the two reservoirs is implemented by the imposition of "Maxwell boundary 
conditions": when a particle hits the left (right) boundary it get reflected such that the velocity 
component parallel to the boundary [v^) derives a new velocity from normal Maxwell distribution 
at the given temperature. The velocity component normal to the boundary {vy) is taken from a 
Maxwell speed distribution given by 

,2 



2 

m 



mv 



corresponding to a temperature '7i('7^) at the left (right) boundary. When the two temperatures 
are different, a net energy flux Je in the y-direction results in the steady state condition, which 
is computed by averaging the kinetic energy added (or removed) by each reservoir per unit time 
and surface. Temperature profiles are obtained from the local kinetic energy density. We use the 
standard velocity Verlet scheme of molecular dynamics with equal time update, except when the 
particles collide with the "hard wall" reservoirs. We treat the collision between the particles and 
the reservoir as that between a hard disk of unit diameter colliding against a hard structure less 
wall. If the time, t^,, of the next collision with any of the two reservoirs at either end is smaller 
than 5t, the usual update time step of the MD simulation, we update the system with t^,. During 
the collision with the walls Maxwell boundary conditions are imposed to simulate the velocity of 
an atom emerging from a reservoir at temperatures Ti or T2. 

The system is allowed to reach the steady state where the current is the same in the entire 
system. To check for local thermal equilibrium , we first define the local temperature T(y) = 
(1/2 mv'^iy)), where the averaging is done locally over strips parallel to the x-direction. This is 
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Figure 6.9: (a) {v'^{y)) and 8T^(?/) as a function of y, the system coordinate perpendicular to 
the reservoirs, for /i = 8.0. (b) T{y) as a function of y, the system coordinate perpendicular to the 
reservoirs, for = 7.0 (Red line ) and n = 8.0 (Blue line). 

called the kinematic temperature. In Fig. 16.9( a) we plot, from the MD simulations with n = 8.0, 
(v^) and see that it overlaps with 8T^ locally. This is consistent with the thermodynamic relation 
{v'^) = 8T^. Thus the system preserves local thermal equilibrium. Note that this verification in 
turn justifies our definition of the local temperature. The temperature profile (Fig. IG.Qf b)) shows 
that (1) the thermal conductivity of the solid region is larger than the liquid, which is expected 
because of the larger density of the former and (2) there is a significant jump in the temperature 
as one crosses the two interface. Such a jump in the temperature is also expected and is due to 



the Kapitza or contact resistance {Rr) 2I1- This is defined as, 

AT 

Rk = — (6.5) 

JE 

where AT is the difference in temperature across the interface. We are next interested in 
determining this Kapitza resistance across the solid liquid interface as a function of the strength 
of the external potential /i. 

We have already seen that with increasing fi, the system shows a jump in the density of the 
solid region corresponding to the addition of an entire layer of atoms. From the profile shown 
in Fig. IG.Qf b). the Kapitza resistance is easily obtained by dividing the temperature jump by the 
energy flux. Note that a slight dependence of on T is visible in Fig. IG.Qf b) with a larger 
temperature jump on the "warm" side. The results shown correspond to the average values of 
Rk between the "warm" and "cold" sides. 

The plot of Rk as a function of fi shows a distinct jump as a layer is included (for a value 
of /i close to 7) in the solid region as shown in Fig. 16.101 The jump in Rk is also accompanied 
by a local peak at the transition. In Fig. 16.11! we have plotted the heat flux through the system 
as a function of /i. As /i increases, the atoms from the surrounding liquid get attracted into the 
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Figure 6.10: Plot of the Kapitza resistance, Rk, as a function of /i, shows a jump at the layering 
transition 



potential well and the density of the liquid progressively becomes lower than that of the solid. 
This results in an overall decrease in the heat flux and consequently the overall conductivity. 
However, close to the layering transition at ~ 7 — 8 there is a local peak in the value of 
the heat flux suggesting that a significant amount of kinetic energy is exchanged between the 
liquid and solid through the interface at the layering transition. The thermal conductivity of the 
solid As also shows a peak at the transition (Fig. I6.12j) . The combined effect of the enhanced 
Kapitza resistance as well as enhanced conductivity of the solid can be summarized by defining 
the Kapitza length Ik, which is the length of solid which has a thermal resistance equivalent to 
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Figure 6.11: Plot of the heat flux, je, as a function of the trap depth, ^u. Note that the overall 
flux decreases as a function of /i. 
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Figure 6.12: Plot of the thermal conductivity of the solid, As as a function of /x. increases as 
the solid takes in a complete layer of atoms. 

the liquid solid interface. This is reminiscent of the "effective liquid" concept which has been 
used in explaining some of the features of acoustic shock absorption in the last section. The 
Kapitza length is given by, 

Ik = RkXs (6.6) 

Note that substitution of the definition of the Kapitza length in the equation for the Kapitza 
resistance (Eqn. 6.3) reproduces Fourier's law once the thermal gradient is identified with AT/Ik- 
The Kapitza length also shows a peak near the layering transition (Fig. 16.13]) . 
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Figure 6.13: Plot of Kapitza length, Ik as a function of fi also shows a jump at the layering 
transition. 
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With the help of these results we may conclude that the layering transition has a profound effect 
on the thermal properties of the liquid solid interface. An important consequence of this study 
is the possibility that the thermal resistance of interfaces may be altered using external potential 
which cause layering transitions in a trapped nano solid. We believe that this phenomenon has 
the potential for useful applications for e.g. as tunable thermal switches or other nano engineered 
devices. 
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